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Preface 


The  main  purpose  of  this  diesis  research  is  to  apply  dithering  techniques  for 
enhanced  idendflability  of  uncertain  parameters  using  moving-bank  Multiple  Model 
Adaptive  Estimation  and  Control  (MMAE/MMAC)  algorithms.  These  techniques  are 
applied  to  the  SPICE  2  flexible  space  structure  in  order  to  quell  induced  structural 
vibrations.  Results  of  this  research  indicate  that  purposefully  constructed  dither  signals 
can  enhance  parameter  identification  to  a  significant  degree. 
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Abstract 

Dithering  techniques  for  enhancing  uncertain  parameter  identification  with  moving- 
bank  Multiple  Model  Adaptive  Estimation  (MMAE)  and  Control  (MMAC)  algorithms  are 
analyzed  in  this  thesis.  The  dithering  techniques  and  multiple-model  adaptive  algorithms 
are  applied  to  the  SPICE  2  flexible  space  structure. 

The  dithering  techniques  studied  include  purposefully  constructed  square  wave, 
sine  wave  (fixed  and  swept  frequency)  and  wide-band  noise  wave  forms.  Purposeful  rigid- 
body  slew  commands  are  also  used  in  order  to  excite  the  structure’s  flexible  bending 
modes.  Dither  inputs  into  the  structure  are  performed  in  an  effort  to  enhance  the  open- 
loop  identifiability  of  the  uncertain  parameter,  namely  a  scalar  multiplier  on  the  undamped 
natural  frequencies  of  the  flexible  bending  modes.  Correct  identification  of  this  parameter 
enables  the  MMAE  and  MMAC  algorithms  to  quell  any  vibrations  induced  into  the 
structure,  in  the  face  of  varying  parameter  values. 

The  results  of  this  study  indicate  that  purposefully  constructed  dither  signals  do 
enhance  parameter  identification  significantly.  The  parameter  identification  enhancement 
due  to  rigid-body  slew  command  input  was  shown  to  be  less  effective  than  that  with 
sinusoid  and  wide  band  noise  dither  inputs,  however. 


MULTIPLE  MODEL  ADAPTIVE  CONTROL  OF  A  LARGE 
FLEXIBLE  SPACE  STRUCTURE  WITH  PURPOSEFUL 
DITHER  FOR  ENHANCED  IDENTIFIABILITY 


1.  Introduction 


Kalman  filters  are  very  useful  in  providing  state  estimates  for  systems  driven  by 
known  inputs  and  white  Gaussian  noise.  In  any  state  estimation  problem,  however,  there 
is  the  underlying  problem  of  obtaining  a  correct  dynamics  model  for  the  filter.  It  may  well 
be  possible  to  obtain  a  model  for  a  typical  system  or  one  for  an  assumed  set  of  conditions, 
but  real  systems  do  not  stay  constant.  The  dynamic  parameters  describing  a  system  will 
undoubtedly  vary  over  time.  (We  term  these  uncertain  parameters).  Tnus,  the  problem  is 
a  single  Kalman  filter  model  may  not  be  robust  enough  to  meet  desired  specifications. 

There  are  many  methods  to  overcome  the  uncertain  parameter  problem,  and  these 
are  thoroughly  discussed  in  [20:  Chapter  10].  One  particular  method  will  be  addressed 
here:  Multiple  Model  Adaptive  Estimation  (MMAE)  [20:  129-136].  With  this  method  are 
developed  various  (multiple)  Kalman  filter  models,  each  based  on  a  different  value  of  a  set 
of  uncertain  parameters.  For  example,  the  uncertain  parameters  for  a  second  order  system 
might  be  the  damping  ratio  and  the  undamped  natural  frequency.  If  these  were  each 
known  to  vary  over  a  particular  range,  then  multiple  filters  could  be  constructed,  each 
with  a  different  value  for  the  damping  ratio  and  undamped  natural  frequency. 
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With  the  multiple  models  developed,  each  is  (conceptually)  used  as  the  basis  for  a 
single  elemental  filter  to  form  a  “bank”  or  column  of  filters.  The  output  of  each  filter  (the 
state  estimate  based  on  a  particular  parameter  value)  is  given  a  probability  that  it  is  the 
estimate  based  on  the  correct  parameter  value.  This  probability  is  calculated  from  a 
history  of  measurement  information  and  from  known  system  characteristics  and  noise 
statistics  (to  which  adaptation  is  not  required).  Each  probabilistically  weighted  estimate  is 
then  summed,  resulting  in  a  probabilistically  weighted  average  of  the  single  filter  state 
estimates. 

The  structure  of  the  Multiple  Model  Adaptive  Estimator  can  be  seen  in  Figure  1-1. 
As  shown,  the  bank  consists  of  K  separate  Kalman  filters,  each  based  on  particular 
parameter  value  from  the  set  {a,,  a2,  .  .  .  ,  a*}.  Each  filter  obtains  the  current 
measurement  input,  Z,  and  from  this  outputs  a  residual,  rt ,  and  a  state  estimate,  xk.  The 
residual  (the  difference  between  the  measurement  and  best  prediction  of  that  measurement 
based  on  previous  data)  is  used  to  compute  the  conditional  probability,  which  is  then  used 
to  form  the  probabilistically  weighted  average  of  the  K  state  estimates. 

The  MMAE  presented  thus  far  is  considered  a  “fixed-bank”  MMAE,  and  as  it 
appears  in  Figure  1-1,  the  MMAE  is  based  upon  a  finite  number  of  parameter  points,  i.e., 
a  finite  number  of  Kalman  filters.  This  finite  size  of  the  parameter  set  used  for  the  MMAE 
is  necessary  because,  if  the  parameter  set  were  to  consist  of  all  possible  points  from  a 
continuous  space,  the  number  of  filters  required  would  be  infinite.  Thus,  only  a  finite 
number  of  points  in  the  appropriate  range  of  values  are  chosen  for  the  parameter  set.  This 
is  then  called  a  discretized  parameter  space.  Now,  even  with  a  finite  number  of  Kalman 
filters  within  this  discretized  parameter  space,  the  number  of  filters  can  be  quite  large  and 
thus  produce  a  considerable  computational  burden. 
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Figure  1  - 1 .  Multiple  Model  Adaptive  Estimator[20: 132] 


To  lessen  this  burden,  the  concept  of  a  subset  of  the  original  parameter  set  is 
introduced.  This  subset  is  a  small  grouping  within  the  parameter  space,  always  centered 
on  the  current  estimate  of  the  true  parameter  value.  This  small  subset  can  then  “move” 
within  the  entire  parameter  space,  depending  on  the  location  of  the  parameter  estimate, 
and  thus  the  term  “moving  bank”  [17].  A  discussion  on  the  moving  bank  will  be  covered 
in  more  detail  in  Section  1.2.3. 
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With  the  state  estimate  in  hand,  we  will  now  discuss  the  controller  to  be  used  m 
this  research.  Corresponding  to  linear  Kalman  filters  driven  by  white  Gaussian  noise,  a 
controller  based  on  a  Linear  model  with  a  Quadratic  cost  criterion  driven  by  white 
Gaussian  noise  (LQG)  will  be  used.  This  stochastic  controller  will  utilize  the  same  K 
state  estimates  produced  by  the  MMAE.  Each  separate  state  estimate  is  multiplied  with  a 
corresponding  controller  gain  of  a  deterministic  optimal  regulator,  -Gc\  to  synthesize  the 

optimal  LQG  controller  for  that  specific  parameter  vector  value.  Just  as  there  are  K 
separate  Kalman  filters,  each  based  on  the  parameter  value,  a*,  there  are  K  separate 
controller  gains,  each  based  on  the  parameter  value,  at.  The  structure  of  this  Multiple 

Model  Adaptive  Controller  (MMAC)  can  be  seen  in  Figure  1-2.  As  shown,  the  controller 
gains  are  cascaded  to  form  control  inputs,  ut .  The  same  probabilities  used  in  the  MMAE 

are  used  here  to  form  a  probabilistically  weighted  average  control  input,  u. 

The  research  involved  in  this  thesis  is  a  direct  continuation  of  the  work  of 
Gustafson  [2,  3],  Gustafson  applied  the  moving-bank  MMAE  and  MMAC  concepts  to  a 
large  flexible  structure  developed  at  Phillips  Laboratory,  Kirtland  AFB,  NM  (see  Figure  1- 
3).  This  structure  is  part  of  a  SPace  Integrated  Controls  Experiment  (SPICE).  Section 
1.2.1  gives  a  brief  description  of  the  SPICE  model.  At  the  time  of  Gustafson’s  work,  the 
second  model  had  been  developed  (SPICE  2).  The  primary  goal  of  this  thesis  is  to 
develop  an  MMAE/MMAC  to  quell  any  vibrations  which  are  induced  into  the  structure, 
despite  uncertain  parameters.  Earlier  contractor  work  [14]  yielded  a  controller  design  that 
would  almost  meet  design  specifications  if  the  parameters  were  assumed  to  be  correctly 
known,  but  which  caused  closed-loop  system  instabilities  if  the  real  parameter  values 
varied  by  as  little  as  2%  from  the  assumed  values.  Gustafson  did  accomplish  this  goal,  but 
the  desired  quelling  specifications  set  by  Phillips  Laboratory  were  not  met.  (Phillips 
Laboratory  requires  the  Line-of-Sight,  or  LOS,  deviation  to  be  within  1  p-radian.)  The 
work  of  this  thesis  will  be  to  apply  enhancements  to  the  existing  models  put  forth  by 
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Gustafson  in  order  to  determine  what,  if  any,  improvement  in  the  MMAE/MMAC  can  be 
made  towards  meeting  the  desired  quelling  specifications.  These  enhancements  include: 

•  Use  dithering  techniques  to  enhance  uncertain  parameter  identifiability,  prior  to 
engaging  closed-loop  control 

•  Investigate  the  effects  of  purposeful  rigid-body  slew  commands  in  regard  to 
enhancing  uncertain  parameter  identifiability 

•  Increase  the  measurement  sensor  precision 

•  Investigate  an  alternative  parameter  space  discretization 
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1.1  Notation 


The  notation  used  in  this  thesis  is  consistent  with  that  of  [19].  Scalars  in  either 
upper  or  lower  case  letters  are  denotes  by  italic  type.  Deterministic  vectors  and  matrices 
are  denoted  by  boldface  roman  lower  and  upper  case  letters,  respectively,  e.g.,  X  and  A. 
Stochastic  vectors  and  matrices  are  denoted  by  boldface  sans  serif  lower  and  upper  case 
letters,  respectively,  e.g.,  X  and  A.  Similar  to  deterministic  vectors,  realizations  of 
random  variables  are  denoted  by  boldface  roman,  e.g.,  X. 


1.2  Background 

This  section  presents  an  introduction  to  four  areas,  the  concepts  for  which  will  be 
fully  developed  in  Chapters  2  and  3.  These  areas  are:  (1)  the  system  model,  (2)  Multiple 
Model  Adaptive  Estimation  (MMAE),  (3)  Moving-Bank  MMAE  and  (4)  Moving-Bank 
Multiple  Model  Adaptive  Control  (MMAC). 


1.2.1  System  Model 

The  full  physical  and  mathematical  description  of  the  SPICE  2  model  will  be 
presented  in  Chapter  3.  Physically,  the  SPICE  2  structure  is  6.19  meters  wide  at  the  base 
and  is  8.14  meters  tall  [14].  In  the  simplest  of  terms,  the  space  structure  can  be  described 
as  a  tripod  structure  with  a  large  mirror  assembly  at  its  base  and  a  small  mirror  assembly  at 
the  top.  In  a  realistic  environment,  the  structure  will  experience  vibrations  (the  three 
tripod  legs  will  vibrate,  causing  the  two  mirrors  to  be  out  of  alignment;  high  frequency 
oscillations  of  the  stiff  mirror  assemblies  will  also  take  place,  but  with  much  smaller 
amplitude  and  effect  on  line-of-sight).  In  order  to  quell  these  vibrations,  actuators  are 
used  which  are  located  on  the  tripod  legs  and  around  the  base  of  the  structure.  In  order 
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for  this  quelling  of  the  structural  vibrations  to  be  accomplished,  the  true  structural 
parameters  must  be  accurately  estimated  (done  by  the  MMAE).  The  parameters  to  be 
estimated  consist  of  the  undamped  natural  frequencies,  ©„,  of  the  structural  bending 

modes  (see  Chapter  2).  Previous  work  using  a  simpler  structural  model  used  the 
undamped  natural  frequencies,  (0n,  and  the  damping  ratios,  £,  of  the  structural  bending 

modes  [1, 2,  3, 5,  7,  8,  1 1, 12,  22,  23,  26, 29,  32].  However,  Gustafson  [2:  Chapter  5  pg. 
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36]  discovered  that,  for  the  SPICE  2  model,  adaptation  to  the  damping  ratio  parameter 
was  not  needed  and  so  is  not  used  in  this  research. 


1.2.2  Multiple  Model  Adaptive  Estimation 

As  indicated  in  the  introductory  paragraphs,  the  need  for  an  adaptive  estimator 
comes  from  the  problem  of  uncertain  parameters  in  the  system  model  from  which  a  single 
Kalman  filter  might  be  formed.  The  multiple  model  approach  utilized  here  involves  a  set 
of  parameters  representing  various  discrete  values  each  parameter  can  take  on  in  a  real 
environment.  For  the  discussion  that  follows,  an  example  parameter  space  is  formed  using 
the  undamped  natural  frequency  and  damping  ratio  as  the  two  uncertain  parameters.  The 
discretized  parameter  space  is  formed  by  assuming  a  range  of  values  for  each  of  the  two 
parameters.  For  example,  let  <Dn  and  £  each  take  on  10  different  values.  Then  100 
combinations  of  ( »„  and  £  would  be  formed  (K  =  100).  In  this  case  there  would  also  be 
100  Kalman  filters,  each  corresponding  to  a  different  combination  of  0>„  and  £.  This 
100-point  parameter  space  can  be  thought  of  existing  in  a  2-dimensional  grid.  This  idea  is 
illustrated  in  Figure  1-4. 

This  two-dimensional  layout  of  the  filters  is  important  because  is  allows  us  to 
visualize  the  K  Kalman  filters  as  a  bank  which  surrounds  the  location  of  the  true 
parameter  value  (or  actually  a  best  estimate  of  that  location).  The  MMAE  algorithm  is 
based  on  utilizing  the  residual  information  from  each  filter  and  the  assumption  that  the 
filter  with  the  smallest  residual  will  most  likely  be  the  filter  with  parameter  values  that  are 
closest  to  the  true  values  [20:  133].  The  residuals,  rk,  are  used  to  produce  a  probability 
weighting  factor.  The  smaller  the  residual,  the  larger  the  probability  which  is  assigned  to 
that  filter’s  output.  Then,  as  indicated  in  the  introduction,  this  probability  is  used  to  form 
the  probabilistically  weighted  state  estimate.  Figure  1-1  illustrates  the  pure  Bayesian 
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Figure  1-4.  Full-Bank  MMAE 
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formulation  of  the  state  estimate.  There  are,  however,  other  methods  of  utilizing  the 
residuals  in  calculating  the  probabilities  and  forming  the  probabilistically  weighted  state 
estimate.  The  four  methods  of  performing  this  are:  (1)  Bayesian,  (2)  Maximum  A 
Posteriori  (MAP),  (3)  Bayesian  with  a  Maximum  Entropy  with  Identity  Covariance 
(ME/I)  computation,  and  (4)  MAP  with  ME/I.  These  four  methods  will  be  fully  discussed 
in  Chapter  2. 


1.2.3  Moving-Bank  MMAE 

The  example  parameter  space  considered  in  this  discussion  requires  100  Kalman 
filters  to  be  run  simultaneously.  Obviously,  this  is  would  be  quite  a  burden  for  any  on-line 
computer  application.  As  indicated  in  the  introduction,  the  concept  of  a  moving-bank  of 
filters  is  introduced  in  order  to  alleviate  this  computational  burden.  The  number  of  filters 
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required  in  this  smaller  bank  is  that  which  will  surround  the  parameter  value.  For  our  two- 
dimensional  parameter  space,  nine  filters  have  typically  been  used.  This  configuration 
allows  a  small  but  reasonable  number  of  filters  to  perform  the  task  of  state  and  parameter 
estimation  with  the  added  benefit  of  having  one  central  point  for  the  basis  of  an  elemental 
filter.  In  general  we  say  the  smaller  moving-bank  consists  of  J  filters  (parameter  points), 
where  J  <  K.  The  arrangement  of  parameter  points  discussed  to  this  point  is  called  finely 
discretized  (all  discrete  points  are  adjacent  to  each  other).  Now,  when  the  true  parameter 
values  vary,  within  the  confines  of  the  defined  parameter  space,  the  small  subset  of  filters 
will  “move”  in  order  to  keep  the  computed  estimate  of  the  true  parameter  value  centered 
within  the  bank.  This  moving-bank  concept  is  illustrated  in  Figure  1-5.  As  the  parameter 
values  change  and  “move”  to  a  new  location  in  the  parameter  space,  so  does  the  bank  of 
filters.  What  is  actually  taking  place  is  the  activating  and  deactivating  the  appropriate 
Kalman  filters,  such  that  the  bank  is  centered  properly. 

The  moving-bank  should  be  able  to  track  a  slowly  varying  parameter  in  this 
manner.  However,  should  the  parameter  “jump”  in  value  a  great  distance  from  its  original 
location  in  the  parameter  space,  a  new  logic  must  be  implemented.  If  the  true  parameter 
value  changes  to  a  location  in  the  parameter  space  outside  the  “bounds”  of  the  finely 
discretized  bank,  all  filters  associated  with  the  finely  discretized  bank  will  have  large 
residuals,  and  the  bank  configuration  must  change  to  a  coarse  discretization,  which  uses 
filters  along  the  outer  edges  of  the  parameter  space.  In  this  way  the  true  parameter  value 
will  be  encompassed  by  the  bank.  Once  the  new  location  of  the  parameter  value  is  found, 
the  bank  can  be  finely  discretized  again.  The  coarse  discretization  concept  is  illustrated  in 
Figure  1-6.  Also  note  that  the  coarse  discretization  can  be  used  during  an  initial 
acquisition  when  the  parameter  value  is  not  known. 
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Figure  1-5.  Moving  Bank  MMAE  Fine  Discretization 


Just  as  there  are  several  methods  of  implementing  the  MMAE,  there  are  several 
methods  for  the  decision-making  logic  of  when  and  in  what  direction  to  move  the  bank, 
and  when  to  expand  and  contract  the  bank.  These  methods  are:  (1)  Residual  Monitoring, 
(2)  Parameter  Position  Estimate  Monitoring,  (3)  Probability  Monitoring,  and  (4) 
Parameter  Estimation  Error  Covariance  Monitoring,  These  methods  will  be  fully 
discussed  in  Chapter  2. 


1.2.4  Moving-Bank  MMAC 

The  LQG  stochastic  controller  will  be  implemented  in  the  form  of  a  regulator  since 
the  desired  quelling  action  necessitates  controlling  all  the  structural  states  to  zero  [18,  21: 
19-20].  Rigid-body  states  are  not  included  in  the  discussion  of  the  Kalman  filter  or  the 
stochastic  controller,  since  the  purpose  of  the  filter/controller  is  to  quell  the  vibrating 
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Figure  1-6.  Moving-Bank  MMAE  Coarse  Discretization 


motion  of  the  structure,  not  the  rigid-body  motion.  A  truth  model  rigid-body  motion  will 
be  discussed  in  Chapters  2, 3  and  4. 

The  LQG  stochastic  controller  used  in  this  research  is  based  on  the  certainty 
equivalence  property  [21:  17].  That  is,  the  controller  gain,  Gc*,  is  formed  assuming  that 
the  optimal  gain  calculated  for  perfectly  known  states  is  the  same  optimal  gain  to  be  used 
with  state  estimates.  Each  control  gain,  as  shown  in  Figure  1-2,  is  calculated  for  each  of 
the  K  parameter  values,  just  as  with  the  K  Kalman  filters  in  the  MMAE  design. 

The  combination  of  the  MMAE  output  with  the  multiple-model  LQG  controller 
can  take  on  many  different  forms.  Figure  1-2  illustrates  a  pure  MMAC  structure  with 
each  Kalman  filter  cascaded  with  a  single  LQG  controller  that  makes  the  same  parameter 
value  assumption  as  the  corresponding  filter.  There  are  slight  variations  of  this  method 
called  “Modified”  MMAC  Control  and  MMAC  with  MAP  control.  The  pure  MMAC  and 
“Modified”  MMAC  methods  are  referred  to  as  Bayesian  because  their  output  is 
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probabilistically  weighted  using  the  hypothesis  probability  calculation,  as  shown  in  Figure 
1-2.  The  MMAC  with  MAP  method  is  a  form  of  MMAC,  but  instead  of  using  each 
conditional  probability  and  elemental  controller  to  form  a  probability-weighted  average  (as 
with  the  Bayesian  form),  this  method  picks  the  probability  with  the  largest  value  and  uses 
the  associated  single  elemental  controller  to  form  the  output  control. 

Other  forms  of  multiple  model  control  exist  and  are  termed  MMAE  based  control. 
These  methods  differ  from  the  MMAC  forms  in  the  way  the  state  estimates  are  utilized 
and/or  in  the  way  the  final  control  output  is  formed.  The  MMAE  based  methods  of 
control  are:  Single  Fixed-Gain  Control,  Single  Changeable-Gain  Control,  and  “Modified” 
Single  Changeable-Gain  Control.  Each  of  these  methods  of  control  will  be  discussed  in 
Chapter  2. 


1.3  Past  Research 

This  section  will  present  a  brief  summary  of  research  accomplished  in  the  areas  of: 
(1)  applying  MMAE/MMAC  methods  to  flexible  structures  and  (2)  applying  dithering 
techniques  to  enhance  parameter  identifiability. 


1.3.1  MMAE/MMAC  Techniques 

Research  was  begun  by  Hentz  [5,  22]  in  the  area  of  multiple  model  filters  applied 
to  the  problem  of  quelling  vibrations  in  a  flexible  structure.  Hentz’ s  research  showed  that 
the  moving-bank,  versus  the  fixed-  or  full-bank,  adaptive  filter  and  controller  were  viable 
solutions.  The  moving-bank  MMAE/MMAC  reduces  the  computational  load  substantially 
while  maintaining  adequate  performance.  Hentz  showed  this  load  reduction  to  be  an  order 
of  magnitude  lower  [5].  Thus,  a  foundation  was  laid  for  continued  research  in  this  area. 
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Filios’  [1]  research  indicated  that  ambiguity  function  analysis  can  provide  insight  as 
to  which  (if  any)  unknown  parameters  need  to  be  estimated.  This  insight  can  be  very 
useful  when  approaching  a  new  system  for  the  first  time. 

Kamick  [7,  8]  applied  the  MMAE/MMAC  methods  to  a  more  complex  model  than 
had  been  done  previously,  namely  that  of  a  two-bay  truss  structure  attached  to  a 
spacecraft.  Kamick’ s  algorithms  had  problems,  though,  correctly  identifying  the  uncertain 
parameters,  due  to  the  sensors  being  too  imprecise.  However,  his  efforts  did  give  accurate 
state  estimation. 

Lashlee  [11,  12]  followed  up  on  the  problems  encountered  by  Kamick.  He  found 
that  with  a  better-tuned  filter  and  more  accurate  measurement  data,  the  parameters  and 
states  could  be  accurately  estimated.  Basically,  if  there  is  too  much  measurement  noise, 
the  multiple  model  algorithm  cannot  distinguish  between  good  and  bad  models  of  the 
uncertain  parameters.  Lashlee  also  found  that  the  parameter  space  discretization  plays  an 
important  role  in  the  parameter  identification  performance. 

Van  Der  Werken  [32]  applied  the  six-state  filter  from  Lashlee’s  research  to  a 
higher-order  “truth  model”  for  simulation  (24-state  truth  model  versus  a  6-state  filter 
model).  In  his  robustness  study,  he  showed  that  the  MMAE  and  MMAC  performed 
poorly  when  comparing  truth  and  filter  model  bending  mode  states  directly. 

Schore  [23,  29]  continued  Van  Der  Werken’s  robustness  study.  His  approach  was 
to  compare  the  truth  and  filter  model  representations  of  the  flexed  space  structure  at 
various  physical  locations  on  the  structure.  That  is,  he  compared  the  position  and  velocity 
(the  shape  of  the  structure)  at  the  two  ends  and  at  the  midpoint  of  the  two-bay  truss,  truth 
versus  filter  [29:  Chapter  1  pg.  22].  Using  this  method,  he  showed  that  the  reduced  order 
filter  model  did  not  seriously  degrade  the  performance  of  the  estimator  or  controller.  In 
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fact,  his  MMAC  maintained  very  adequate  control  in  the  face  of  an  unmodeled 
disturbance. 

Moyle’s  [26]  research  concentrated  on  streamlining  the  software  that  had  been 
developed  thus  far  and  on  further  filter  and  controller  tuning  to  enhance  performance. 
Moyle’s  research  showed  marked  performance  improvement  of  the  MMAE/MMAC,  when 
compared  against  an  artificially  informed  benchmark  (an  algorithm  that  was  artificially 
provided  knowledge  of  the  true  parameter).  This  performance  was  obtained  by 
implementing  a  Maximum  Entropy  with  Identity  (ME/I)  covariance  computation  approach 
(see  Section  2.3.2)  to  parameter  estimation  along  with  the  “Modified”  MMAC  algorithm 
for  which  the  filter  and  controller  were  well  tuned  [26:  Chapter  1  pp  13-15].  The  ME/I 
approach  and  the  “Modified”  MMAC  algorithm  will  be  fully  discussed  in  Chapter  2. 

The  most  recent  work  is  that  of  Gustafson  [2,  3],  who  applied  the  current  research 
to  the  SPICE  2  model.  This  was  the  first  time  this  line  of  research  was  applied  to  a  total, 
actual  system.  Gustafson’s  work  showed  good  performance  of  the  moving-bank 
MMAE/MMAC  with  a  134-state  filter  model  versus  a  194-state  truth  model.  However, 
with  more  stringently  reduced  order  filter  models,  the  MMAC  could  not  quell  the 
disturbance  vibrations  adequately.  Additionally,  the  good  performance  of  the  controller 
was  marred  by  its  inability  to  keep  the  LOS  deviations  within  the  1  (i-radian  specification. 
The  inability  to  meet  this  specification  was  attributed  to  poor  actuator  and  sensor 
performance  and  did  not  lie  in  the  MMAE/MMAC  adaptation  process:  even  a  nonadaptive 
artificially  informed  controller  failed  to  meet  the  specification. 

Gustafson  also  ran  tests  using  an  expanded  H  matrix  which  resulted  in  better 
quelling,  but  the  specification  was  still  not  met  [2:  Chapter  5  pg.  3].  This  expanded  H 
matrix  simulated  relative  position  and  velocity  measurements  of  the  sensor  measurement 
points  on  the  structure.  These  extra  measurements  were  not  physically  realizable  but  were 
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useful  for  a  comparison  study.  The  performance  improvement  shown  here  can  be  equated 
to  an  expected  equal  or  greater  performance  gain  if  the  sensor  measurement  precision  is 
improved. 


1.3.2  Dithering  Applications 

Much  work  has  been  done  in  the  past  using  dither  signals  to  enhance  parameter 
identifiability.  A  few  of  the  applications  which  utilized  these  techniques  will  be  presented 
here. 

In  this  direct  line  of  research  prior  to  Gustafson  [1,  5-8,  11,  12,  22,  23,  26,  29, 
32],  dithering  signals  were  input  to  the  system  in  order  to  obtain  an  adequate  uncertain 
parameter  identification;  there  were  no  specifically  designed  dithers,  just  some  rather 
arbitrarily  chosen  excitation  signals.  This  was  necessary  because  the  dynamics  driving 
noise  was  too  small  for  the  parameters  to  be  identified  without  some  form  of  forced 
excitation.  The  system  needed  to  be  “shaken  up”  with  a  dither  signal  in  order  for  any 
parameter  identification  to  take  place.  TypicalK  square  wave  dither  signals  were  used  in 
this  previous  work,  but  with  no  purposeful  structure.  However,  with  the  SPICE  2  model 
used  by  Gustafson,  a  sufficient  magnitude  of  dynamics  driving  noise  was  present  to 
operate  somewhat  successfully  in  open  loop  without  the  use  of  a  dither  signal,  so  none 
was  applied.  Gustafson  showed  that  the  MMAE  was  somewhat  able  to  track  a  varying 
parameter  under  these  conditions  [2:  Chapter  5  pp.  46-49]. 

Menke  and  Maybeck  [24,  25]  used  a  variety  of  dither  input  signals  in  an  MMAC 
application  for  detecting  sensor  and/or  actuator  failures  in  a  flight  control  system:  sine 
wave,  square  wave,  triangle  wave,  and  pulse  trains,  as  well  as  purposeful  maneuver 
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commands.  All  of  these  dither  inputs  performed  well,  but  the  sine  wave  gave  the  best 
performance  [24:  203]. 

Hanlon  [4:  36-37]  looked  at  the  eigenvalues  in  a  flight  control  system  to  determine 
the  best  sine  wave  dither  signal  to  use  for  enhancing  failure  detection.  The  eigenvalues 
provided  direct  insight  as  to  which  dither  frequency  to  use.  Before  using  the  insight 
gained  from  the  eigenvalues,  an  arbitrarily  chosen  dither  frequency  was  chosen  which  did 
not  produce  useful  results. 


1.4  Objectives 

The  purpose  of  this  thesis  is  threefold.  First,  various  dithering  techniques  will  be 
implemented  in  order  to  enhance  uncertain  parameter  identifiability  for  the  SPICE  2  model 
in  open-loop  estimation  prior  to  closed-loop  control.  It  is  envisioned  that  this  will  provide 
enhanced  initial  parameter  estimates  when  the  closed-loop  control  system  is  turned  on; 
there  is  enough  excitation  in  the  closed-loop  phase  that  the  dither  would  only  be  used  in 
the  prior  open-loop  phase.  The  research  will  also  explore  whether  or  not  naturally 
occurring  slew  commands  will  excite  the  system  sufficiently  without  dithering.  Second, 
increased  sensor  measurement  precision  will  be  investigated  in  an  effort  to  increase 
controller  performance  (decrease  LOS  deviations).  It  is  anticipated  that  reduced-order 
design  models  will  be  more  fruitful  once  measurement  precision  is  improved.  Third,  an 
alternative  to  Gustafson’s  discretized  parameter  space  will  be  developed  in  an  order  to 
improve  MMAE/MMAC  overall  performance. 
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1.5  Scope 


Since  this  thesis  research  is  a  direct  continuation  of  Gustafson’s  [2],  this  research 
will  utilize  the  existing  SPICE  2  model  in  addition  to  the  Kalman  filter  and  controller 
tuning  developed  by  Gustafson,  except  where  modifications  are  necessary  (see  next 
section).  The  single  unknown  parameter  (undamped  natural  frequency  of  the  structural 
bending  modes)  will  continued  to  be  used  in  this  research.  The  same  discretized 
parameter  space  developed  by  Gustafson  will  also  be  used,  except  when  a  new 
discretization  is  formed. 

The  limits  placed  on  the  modeling  are  as  follows.  The  truth  model  is  a  linearized 
model  for  the  SPICE  2  structure.  The  Kalman  filter  is  based  on  a  reduced-order  linear 
time-invariant  model,  driven  by  stationary  noise.  The  LQG  controller  is  developed  from 
constant  state  and  control  weighting  matrices  with  no  cross-weighting  matrix.  Gain 
transients  in  both  the  filter  and  controller  will  be  ignored,  in  order  to  exploit  the  steady- 
state  constant-gain  Kalman  filter  in  the  MMAE,  and  the  steady-state  constant-gain  LQG 
controller  in  the  MMAC. 


1.6  Approach 

To  apply  the  dithering  techniques,  two  items  will  be  accomplished.  First,  several 
forms  of  dither  signal  inputs  will  be  considered  in  order  to  observe  the  effects  each  has  on 
enhancing  parameter  identification.  The  different  inputs  to  be  considered  are:  (1)  square 
wave,  (2)  sine  wave,  (3)  swept  frequency  sine  wave,  (4)  wide-band  noise  (limited  to  the 
expected  band  of  undamped  natural  frequencies).  Second,  rigid-body  states  will  be 
augmented  to  the  SPICE  2  structural  model  to  simulate  a  commanded  rotation  or  slewing 
motion  of  the  structure.  The  effects  of  this  rigid-body  motion  will  be  observed  with 
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respect  to  how  it  affects  parameter  identification  through  excitation  of  the  bending  modes, 
just  as  with  the  dither  signal  inputs.  Third,  simulations  will  be  run  with  the  sensor 
measurement  precision  increased  by  one,  two  or  more  orders  of  magnitude.  With  an 
increase  in  the  sensor  measurement  precision  (decrease  in  sensor  noise),  the  effects  in  the 
MMAC  performance  will  be  observed  when  based  on  more  stringently  reduced-order 
models  (with  associated  decrease  of  computational  loading). 

With  the  dithering  inputs,  a  comparison  will  be  accomplished  among  the  various 
techniques  to  determine  which,  if  any,  is  “best”.  Minor  modifications  to  the  existing 
software  [2]  will  be  required. 

Increasing  the  sensor  measurement  precision  will  require  retuning  of  the  existing 
Kalman  filter  designs  (assumed  dynamics  noise  strength,  Q,  and  assumed  measurement 
noise  covariance,  R ).  Before  changing  the  precision  by  one  or  two  orders  of  magnitude, 
the  effects  of  perfect  measurements  (no  measurement  noise)  will  be  determined.  In  this 
way,  a  limiting  case  in  the  MMAC  performance  will  be  known,  with  regard  to  increased 
sensor  measurement  precision. 

The  change  to  the  parameter  space  discretization  will  be  either  to  change  the 
location  of  discrete  point  values,  or  to  allow  a  different  number  of  points,  or  both.  The 
existing  span  of  values  for  the  undamped  natural  frequencies  will  be  maintained,  however, 
for  the  purpose  of  making  an  accurate  comparison  of  the  MMAE/MMAC  performance 
with  the  old  and  new  parameter  space  discretization.  Minor  software  modifications  to  the 
existing  software  will  be  required  here  as  well. 
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1.7  Summary 

A  general  introduction  and  background  in  the  area  of  utilizing  the  Kalman  filter 
and  LQG  controller  with  a  multiple  model  approach  to  achieving  desired  performance 
despite  parameter  uncertainty  were  presented  in  this  chapter.  The  techniques  to  be  applied 
toward  this  research  were  also  presented,  namely:  dithering  techniques  for  uncertain 
parameter  open-loop  identification  enhancement,  increased  sensor  measurement  precision 
for  improved  MMAE/MMAC  performance,  and  modified  parameter  discretization  for 
better  parameter  tracking  performance. 

To  motivate  the  desire  for  the  use  of  the  multiple  model  approach  and  the  dithering 
techniques  to  be  implemented,  a  summary  of  previous  research  in  these  areas  was  also 
presented.  The  specific  objectives  of  this  research  along  with  a  defined  scope  and  general 
approach  (methodology)  were  then  presented. 

Chapter  2  will  develop  the  fundamental  concepts  of  the  Kalman  filter  and  LQG 
controller.  These  concepts  will  then  be  used  to  develop  the  moving-bank  MMAE  and 
MMAC  algorithms.  Chapter  2  will  also  present  a  section  on  applying  dithering  input 
signals  and  a  section  on  mathematical  modeling  of  the  physical  structure.  Chapter  3  will 
present  a  summary  of  the  SPICE  2  structure  model  developed  by  Gustafson  [2:  Chapter 
3].  Additionally,  Chapter  3  will  include  the  development  for  augmenting  rigid-body  states 
to  the  SPICE  2  model.  The  software  requirements  and  simulations  to  be  performed  will 
be  presented  in  Chapter  4.  Chapter  5  will  include  a  presentation  of  the  results  found  in 
this  research.  Chapter  6  will  finish  with  a  discussion  of  the  conclusions  drawn  from  the 
results,  along  with  recommendations  for  future  research. 
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2.  Development  of  Fundamental  Concepts 


2.1  Introduction 

Three  main  concepts  will  be  developed  in  this  chapter:  the  Kalman  filter,  the  state 
and  parameter  estimator,  and  the  stochastic  controller.  The  Kalman  filter  development 
will  be  thorough  but  basic.  For  the  state  and  parameter  estimator,  Multiple  Model 
Adaptive  Estimation  (MMAE)  algorithms  will  be  implemented,  and  this  chapter  will 
present  different  methods  of  implementing  the  MMAE.  The  stochastic  controller  will 
utilize  Linear-Quadratic-Gaussian  (LQG)  techniques  as  part  of  the  Multiple  Model 
Adaptive  Controller  (MMAC).  As  with  the  MMAE,  various  methods  of  implementing  the 
MMAC  will  be  addressed. 

In  addition  to  these  three  main  concepts,  two  additional  items  will  be  discussed. 
First,  algorithms  for  implementing  the  dithering  techniques  for  enhanced  parameter 
identification  will  be  developed.  Second,  a  mathematical  modeling  section  will  be 
presented  in  which  physical  and  modal  forms  are  developed  in  order  to  demonstrate  the 
modal  reduction  technique  of  producing  a  reduced-order  filter  model  from  a  truth  model. 

The  development  presented  in  this  chapter  is  independent  of  the  particular  system 
of  concern,  e.g.,  a  flexible  space-structure.  For  this  reason,  the  equations  and  algorithms 
developed  here  can  be  applied  to  almost  any  situation  which  is  constrained  by  the  same 
assumptions  contained  in  this  work. 

2.2  Kalman  Filter  (Bayesian  Estimation) 

In  order  to  utilize  the  simplest  form  of  the  Kalman  filter,  the  system  under  study  is 
assumed  adequately  modeled  by  a  continuous-time,  stochastic  linear  differential  equation 
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driven  by  zero- mean  white  Gaussian  noise.  The  sampled-data  measurement  model  has 
similar  assumptions,  i.e.,  linear  measurements  disturbed  by  zero-mean  white  Gaussian 
noise.  Thus,  the  Kalman  filter  is  based  on  the  following  linear  stochastic  differential 
equation: 


x(f )  =  F(f)x(f ) + B(f  )u(r) + G(f)w(f)  (2.1) 

where  X(f)  is  the  n-dimensional  state  vector,  u(f)  is  the  r-dimensional  deterministic 
control  input  vector,  and  W (f)  is  the  5-dimensional  white  Gaussian  noise  vector.  F(f)  is 
the  n-by-n  dimensional  system  dynamics  matrix,  B(f)  is  the  n-by-r  dimensional 
deterministic  input  matrix,  and  G(f)  is  the  n-by-s  dimensional  noise  input  matrix. 

As  indicated  above,  the  driving  noise  has  statistics  of: 

E {w(f)}  =  0  (2.2) 

f{w(r)w(r'  )T}  =  Q(r)8(f  - 1 )  (2.3) 

where  Q(f)  is  an  s-by-s  dimensional,  symmetric,  positive  semi-definite  matrix  and 
indicates  the  strength  of  the  dynamics  driving  noise;  8(f)  is  the  Dirac  delta  function. 

As  with  all  stochastic  linear  differential  equations,  initial  conditions  must  be 
specified.  These  are  in  the  form  of  x0  and  P0.  With  these  two  initial  conditions  known, 

the  entire  corresponding  Gaussian  probability  density  function  can  be  described.  The 
initial  conditions  are  as  follows: 


£{x(f0)}  =  x0 


(2.4) 


£{[x(f0)-x0][x(f0)-x0]T}  =  P0 


(2.5) 
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Here,  x0  is  the  mean  or  “best  guess”  of  the  initial  state  vector.  P0  is  the  initial  covariance 
or  indicator  of  the  “spread”  of  the  density  function. 

With  the  ultimate  implementation  of  the  Kalman  filter  on  a  digital  computer,  we 
need  a  discretized  set  of  Kalman  filter  equations.  It  is  well  known  that,  in  general,  it  is 
preferable  to  develop  a  discrete-time  filter  around  an  equivalent  discrete-time  system 
model,  rather  than  discretizing  a  continuous-time  filter  [19:261].  Thus,  the  discretized 
system  model  is  given  as: 


x(r, )  =  0(r, ,  )x(rI_1 )  +  Bd  )u(rM )  +  Gd  )wd  )  (2.6) 

where  Qih’h-i)  is  the  state  transition  matrix  associated  with  F(r),  i.e.,  the  solution  to: 

O(r,rM)  =  F(0<l>(/,/1.1)  (2.7) 

with  the  initial  condition: 

0(/,.1,r,_1)  =  I  (2.8) 

Assuming  time  invariant  F,  the  Laplace  transform  solution  yields 

)  =  <&«,  -  tM)  =  =  l  -1  {[SI  -  F]-‘  ^  (2.9) 


If  we  then  assume  constant  u  over  each  sample  period,  Bd  )  is  solved  by: 


B  *(Im) 


fb(^,T)B(x)<fx 
A- 1 


(2.10) 
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Similar  to  the  continuous  time  case,  the  statistics  of  the  dynamics  driving  noise  are  given 
by: 


£{Wd(',-i)}  =  0  (2.11) 

£{*„(', =  Q„(rw)8(WM  (2.12) 

where  5 (i_1)7  is  the  Kronecker  delta  function.  Qd(f)  is  given  by: 

Qd(^-i)=  f  »  t)G(t)Q(t)Gt  (t)Ot  (r, ,  x)dx  (2.13) 

Measurements  are  obtained  by  the  discrete-time  measurement  model: 


Z(*, )  =  H(f,  )X(f, )  +  V(f, )  (2.14) 

where  Z(t, )  is  an  m-dimensional  discrete-time  measurement  vector.  H(ff),  is  an  m-by-n 
dimensional  measurement  matrix  and  V(f(),  is  an  m-dimensional  noise  vector  representing 

the  uncertainty  in  the  measurement  information.  As  with  the  dynamics  driving  noise, 
Wd(fM)  in  Equation  (2.6),  this  vector  is  modeled  as  white,  Gaussian  noise  with  statistics 

given  by: 


£{v(O}  =  0 
^{v(t,)v(r;)T}  =  R(f,)5| 


(2.15) 

(2.16) 


where  R(f, )  is  an  m-by-m  dimensional,  positive  definite,  symmetric  matrix  and  by  is  the 
Kronecker  delta  function.  R(f,  )  being  positive  definite  indicates  that  each  measurement 
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component  is  noise-corrupted  and  that  no  linear  combination  of  these  components  is 
noise-free. 

One  final  assumption  is  that  X(f0),  W(/)and  V(f,)  are  assumed  independent  and 
thus  uncorrelated  with  each  other.  The  equivalence  of  independence  and  uncorrelatedness 
comes  from  the  joint  Gaussian  assumption  [19:205]. 

Now  that  the  discrete-time  system  is  described,  the  propagation/update  cycle  of 
the  Kalman  filter  will  be  addressed.  The  propagation  portion  of  this  cycle  consists  of 
propagating  (moving  ahead  in  time)  the  most  recent  state  estimate  information  at  time, 
t*_  j,  to  time,  t~.  The  superscript  “+”  denotes  the  time  just  after  a  measurement  update, 

while  the  superscript  denotes  the  time  just  prior  to  a  measurement  update.  Just  after 
the  propagation  comes  a  measurement  update.  Thus,  this  cycle  can  be  described  by  the 
following  sequence: 


»*2  * . >*i  »•••*! 


/+!**•• 


The  propagation  equations  are  given  by: 


x(f,~ )  =  )x(£, )  +  Bd  (fw  )u  (fw )  (2. 17) 

P(^,~ )  =  *  h-\  )P(Ci  )^T  (f« » h-\ )  +  Gd  (fw  )Qd  (fM  )GdT  (*,_, )  (2.18) 

The  following  set  of  update  equations  are  used  to  update  the  state  estimate  with  the 
information  provided  by  the  current  measurement: 

K(i, )  =  P«,- )Ht(i,  )[H(i,  )P(tr  )Ht((,  ) + R(r, )]~'  (2.19) 

id* )  =  xdr )  ■ +  K«,  )[z,  -  H(r,  )]  (2.20) 

Pd*)= P dr )  -  K«,  )H«,  )P(r,- )  (2.21) 
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where  z,  is  the  measurement  to  be  incorporated  at  time  f,  (a  specific  realization  of  z(/,)), 
and  K(/, )  is  the  Kalman  filter  gain  matrix.  The  residual  equation  is  given  by: 

r(0  =  zi-H(ri)x(/-)  (2.22) 

where  r(ti )  is  termed  the  residual  and  is  simply  the  bracketed  quantity  in  Equation  (2.20) 
[10:217-218].  The  residual  term  is  singled  out  because  of  its  special  importance  in  the 
system  analysis  and  the  MMAE  to  follow  in  the  next  section.  The  residual  is  seen  to  be 
the  difference  between  the  current  measurement  and  the  filter's  best  “guess”  of  what  that 
measurement  should  be.  This  difference  “guides”  the  filter  in  the  correct  direction  (+  or  -) 
for  updating  the  state  estimate.  The  magnitude  of  the  update  is  controlled  by  the  gain 
matrix.  The  filter  anticipates  that  the  residual  will  be  Gaussian,  of  mean  zero  and  of 
covariance 


A(«i )  =  H«,  )P«  -  )Ht«,)  +  R(!,)  (2.23) 

which  is  inherently  computed  by  the  filter  in  the  gain  calculation  of  Equation  (2. 19).  Since 
x(r(+)  is  formed  from  x(t~)  in  Equation  (2.20),  the  state  estimate  always  incorporates 

knowledge  (information)  of  previous  measurements.  This  is  the  underlying  idea  behind 
recursive  Bayesian  estimation. 

With  the  Kalman  filter  system  of  equations  presented,  we  will  discuss  some  further 
simplifying  assumptions  which  will  minimize  the  computational  work  involved  with 
implementing  such  a  system.  The  following  simplifications  revolve  around  the  current 
assumptions  of  a  time-invariant  system  with  stationary  noise  processes.  With  this  type  of 
system  the  filter  is  known  to  have  a  transient  followed  by  a  steady  state  response  in  P  and 
K  [19:224].  If  the  transient  period  is  relatively  short,  relative  to  the  total  time  of  interest, 
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then  the  steady  state  values  may  be  used  (if  performance  degradation  is  sufficiently  small). 
Thus,  with  the  time-invariant  assumption  and  a  fixed  sample  period,  P“,  P+  and  K  may 
be  precomputed  once  for  all  time,  as  can  Bd,  and  relieve  a  considerable  computational 
burden.  With  the  addition  of  the  stationary  noise  assumption,  Qd  can  be  precomputed 
once  for  all  time  as  well.  For  the  system  under  consideration,  another  valid  assumption  is 
for  R  to  be  constant  as  well. 

With  the  simplifying  assumptions  in  place,  the  last  step  in  the  Kalman  filter  system 
development  is  that  of  “tuning”.  This  process  involves  adjusting  or  “tuning”  the  values  for 
x0,  R  and  Qd .  The  development  of  these  and  the  other  system  model  matrices  (choice  of 
state  variables  in  X;  F,  O,  Bd,  u,  Gd  and  H)  will  be  presented  in  Chapter  3.  The  actual 
tuning  process  will  be  presented  in  Chapter  4. 

2.3  MMAE 

The  basic  understanding  of  the  Multiple  Model  Adaptive  Estimator  (MMAE)  was 
discussed  in  Chapter  1.  This  section  will  expand  on  this  idea,  starting  with  the  uncertain 
parameter  problem  and  the  Bayesian  formulation  of  the  MMAE.  This  is  followed  by  a 
discussion  on  underlying  problems  with  the  basic  or  raw  form  of  the  Bayesian  MMAE. 
This  section  is  concluded  with  a  discussion  on  proposed  enhancements  to  address  these 
problems. 

2.3.1  Bayesian  MMAE 

Uncertainties  are  already  a  given  with  a  stochastic  system  (dynamics  and 
measurement  noise),  but  these  noises  are  white  stationary  processes.  The  unknown  or 
uncertain  parameter  comes  into  play  when  forming  a  system  matrix  (O  in  our  case). 
Note,  that  in  general,  uncertain  parameters  can  affect  Bd,H,Qd  and/or  R  as  well.  The 

Kalman  filter  for  a  particular  system  is  only  as  accurate  as  the  system  matrices  which  make 
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up  the  filter.  If  there  exist  some  uncertainty  in  the  system  at  hand  (damping  ratio  and 
undamped  natural  frequency  of  various  structural  bending  modes  in  our  case),  one  method 
of  handling  this  is  the  MMAE. 

To  understand  the  Bayesian  formulation  of  the  MMAE,  consider  the  following. 
Let  a  be  the  random  vector  composed  of  the  uncertain  parameters.  In  this  initial 
discussion,  a  will  be  considered  a  time-invariant  random  vector.  As  presented  in  Chapter 

1,  the  parameter  space  will  be  comprised  of  a  set  of  K  parameter  vector  values, 
{a,,a2,...at,...aA:}.  This  discretized  space  is  preferred  rather  than  a  continuous  space  for 

obvious  dimensionality  problems;  later  it  will  be  seen  to  allow  an  MMAE  implementation 
with  a  finite  number,  rather  than  an  infinite  number,  of  elemental  Kalman  filters.  (The 
actual  parameter  space  discretization  will  be  developed  in  Chapter  3.)  Now  consider  a 
joint  density  function  of  X(/,)  and  a,  conditioned  on  measurements  Z (f,).  Here,  Z(ff)  is 
the  measurement  history  through  time  f, ,  made  up  of  the  components 
Z(f1),  Z(r2),  Z(r,).  Maybeck  [20:76,129)  states  this  is  the  logical  choice  for  a 

Bayesian  formulation  and  is  given  by: 


(2.24) 


Note  that  Bayes'  rule  [19:81]  is  applied  to  obtain  the  equality.  TUs  form  of  the  density 
function  is  useful  since  the  two  terms  on  the  right  of  Equation  (2.24)  are  readily  known. 
The  first  term  on  the  right  is  the  Gaussian  density  defined  by  the  Kalman  filter  based  on 
the  assumption  that  a  assumes  the  value  a,  with  mean  x(f,+  )  and  covariance  P(f,.+  ) 

[20:129].  The  second  density  function  in  Equation  (2.24)  is  formed  by  considering  the 
probability  of  a  taking  on  the  value  a  at  time  f, .  The  density  function  is  then  described 

by  [20:130]: 
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(2.25) 


K 

falZOt)  (al  ^1 )  =  ^  Pk  (*i  )$(a~ak) 

*= 1 

Then,  for  any  time  f(  ,  the  conditional  hypothesis  probability  is  recursively  defined  as 
[20:131]: 


Pk  (h ) s  Prob{a  =  a*  I2(r,. )  =  Z, } 


fiiUjaZiU-O  (z,  |  at ,  ZM)  pk 

^/z(ti)la^(f._,)(z.- 1  a/»  Z,_i)  Pj  (/M) 
>1 


(2.26) 


where  each  measurement  conditional  density  is  evaluated  as  [20: 132]: 


/ l(rj)|«,Z(»,.1)(Z<  a* » ^i-1 )  = 


(27t)^|At(rl.)| 


1/2 


«p[-Tr.T(',)A;,((,)r,(/,)]  (2.27) 


where  rk(ff)  and  A t(ff)  are  computed  in  the  itth  Kalman  filter  (based  on  the  assumption 
that  a  =  at)  as: 


and 


rt«,)  =  z<f,)-H, («,)*,(/,-) 

(2.28) 

A,«, )  =  H,  (t,  )Pt  («,. ) + Rt  (4  ) 

(2.29) 
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From  this  Bayesian  formulation  we  will  have  a  state  estimate  (conditional  mean)  and 
covariance,  both  of  which  are  conditioned  on  the  past  measurement  history.  The 
conditional  mean  is  formed  by  [20: 131]: 


x(r(+)  =  £{x(ri)IZ(0  =  Z1.} 

=  JL,  /x<J,)l*LZ(r,)(^la*  *  ^Pk  (*«  ) 


L  t=l 


4 


(2.30) 


*=i 


where  xt(f(+)  is  produced  by  the  Mi  Kalman  filter  based  on  the  parameter  at.  Note  the 
difference  between  x(t,+  )  and  xt(f,+).  x(f(+)  is  merely  the  probabilistically  weighted 
average  of  each  of  the  Kalman  filter  state  estimates.  Refer  back  to  Figure  1-1  for  a 
graphical  depiction  of  this  process,  also  known  as  the  multiple  model  filtering  process 
[20:131]. 

If  it  is  desired  to  obtain  the  conditional  covariance  of  X(f, ),  it  can  be  found  from 
[20:131]: 


p«,* )  =  e{[x((,  )-  i(r,*  )Jx(i, )  -  x(f,*  )]T|  Z«, )  =  Z,  ■ 

= i>, <«, ){p, (',*)+[*. (o - *(ol*,  ('/ ) - *(«/)]' }  <2-3i> 
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where  P \(t*)  is  the  covariance  of  each  Kalman  filter’s  state  estimate  based  on  the 
parameter  at,  as  explicitly  precomputable  in  the  fcth  Kalman  filter.  Note  here  that  P(/,+  ) 
cannot  be  precomputed  since  it  is  directly  tied  to  knowledge  of  the  measurement  history. 

In  order  to  evaluate  the  performance  of  the  multiple  model  parameter  estimation,  it 
will  be  useful  to  calculate  the  conditional  mean  and  covariance  of  the  parameter.  These 
are  given  by  [20:132-133]: 


and: 


a(t,)  =  £{a(ri)lZ(r,)  =  Z,}  =  jaf^iwZ^da 

=  JL  “Eft  ('/  )5(“  “  a* ) 

K 

=S**p*(A) 

*=1 

Pi(»,)=£{[a-i((,)Ia-a(0]TIZ(f,)=Z,j 


(2.32) 


(2.33) 


*=! 


2.3.2  Enhancements  to  Bayesian  Estimation 

The  preceding  section  gave  the  Bayesian  formulation  in  multiple  model  form  for 
the  problem  of  state  estimation  (and  parameter  estimation)  in  the  face  of  uncertain 
parameters.  This  section  will  discuss  alternate  formulations  of  the  multiple  model’s  state 
estimate  and  other  enhancements  to  include  in  the  multiple  model  algorithm. 

The  Bayesian  setup  of  the  multiple  model  adaptive  estimation  algorithm  in  its  pure 
form  has  two  potential  drawbacks.  For  each,  there  exist  modifications  or  enhancements  to 
overcome  these  drawbacks.  In  order  to  discuss  the  first  drawback,  consider  the 
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conditional  probability  (tensity  function  in  Equation  (2.27).  Normally  the  density  function 
with  the  largest  magnitude  is  the  one  which  contains  the  smallest  residual  term,  r* .  This  in 

turn  produces  the  highest  probability  for  the  kth  Kalman  filter  in  Equation  (2.26).  This  is 
how  the  algorithm  should  work,  since  the  filter  with  the  smallest  residual  is  likely  to  be 
based  on  a  model  that  is  closest  to  the  truth  model.  However,  if  the  situation  arises  in 
which  each  r/A-lrt  is  of  the  same  relative  magnitude  (which  would  be  desired  to  yield 

equal  probabilities  for  each  elemental  Kalman  filter),  then  the  magnitude  of  the  conditional 


density  function  will  be  driven  by  the  r-^-r  term  in  Equation  (2.27).  In  this  case  the  filter 


with  the  smallest  |At|  will  produce  the  largest  density  value  and  thus  the  largest 
probability.  This  is  certainly  not  the  criterion  upon  which  to  base  the  probabilistically 
weighted  state  estimate.  To  overcome  this  possible  shortcoming  of  the  pure  Bayesian 
formulation,  an  alternative  or  enhancement  is  available  and  is  discussed  below. 

For  this  enhancement  we  will  assume  that  the  residuals  are  Gaussian  with  an 
identity  covariance  matrix.  This  implies  that  the  residuals  will  follow  a  “maximally 
noncommittal  residual  distribution"  [10:26].  The  enhancement  is  termed  Maximum 
Entropy  with  Identity  Covariance  (ME/I).  With  this  method  the  Ak  term  in  the 

conditional  density  function  of  Equation  (2.27)  is  replaced  with  the  identity  matrix,  I . 
The  resulting  density  function  is  given  by: 


la* .  Z,_, )  =  exp[-  ‘  rk  (r,  )rk  (t, )]  (2.34) 

Note  that,  without  this  modification,  we  assume  the  residuals  are  Gaussian  with 
covariance  Ak  and  that  “a  Gaussian  distribution  is  the  maximum  entropy  distribution” 

[29:32,  10:24].  Thus,  the  method  of  using  the  conditional  density  function  with 
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covariance  Ak  (Equation  2.27)  is  termed  Maximum  Entropy  with  Covariance  Ak 
(ME/A). 

A  second  problem  can  occur  if  a  probability,  pk  (ti ) ,  is  ever  calculated  to  be  zero. 
Because  of  the  iterative  nature  of  the  probability  calculation  in  Equation  (2.26),  if  pk  (ti ) 
ever  becomes  zero  it  will  remain  at  zero  for  all  time  thereafter.  This  condition  is  often 
termed  “lock-out.”  To  prevent  such  a  condition,  an  ad  hoc  technique  of  “lower  bounding” 
pk(ti)  is  used  [20:135,  17].  This  enhancement  of  setting  a  lower  bound,  pmn,  on  pk( f.) 

prevents  the  probability  from  ever  becoming  zero.  The  lower  bound  is  set  sufficiently  low 
so  as  to  not  bias  the  computation,  and  at  each  sample  time,  ti ,  after  each  probability  is 

K 

computed  with  Equation  (2.26),  each  pk  (t, )  is  rescaled  such  that  ^pt(/,)  =  1. 

t=i 

In  addition  to  the  enhancements  presented  above,  there  is  another  method  of 
computing  the  state  and  parameter  estimate.  This  method  is  termed  Maximum  A 
Posteriori  (MAP)  and  assumes  the  filter  with  the  highest  associated  probability  will 
provide  the  best  state  estimate.  This  method  can  be  used  with  or  without  the  ME/I 
technique  of  computing  the  probability  density  function. 

2.4  Moving-Bank  MMAE  Development 

The  moving-bank  concept  was  introduced  in  Chapter  1  and  is  an  enhancement  to 
fixed-bank  (full-bank)  MMAE.  In  order  for  the  MMAE  to  be  robust,  the  parameter  space 
must  be  “large”,  but  this,  of  course,  can  create  a  burden  for  the  microprocessor  if  the 
discrete  point  values  are  finely  distributed  throughout  the  space.  To  alleviate  this  burden, 
a  subset  of  J  filters  is  used  at  any  given  time  (with  J  <K).  Then,  for  the  MMAE  to  be 
effective,  the  subset  of  filters  must  “move”  through  the  parameter  space  (to  keep  up  with 
changes  in  true  parameter  value,  a).  In  order  for  the  move  to  work  properly,  the  bank 
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will  also  have  to  be  expanded  and  contracted  periodically.  Additionally,  as  the  bank 
moves  through  the  parameter  space,  new  elemental  Kalman  filters  will  have  to  be 
initialized.  In  this  section,  techniques  for  moving,  expanding  and  contracting  the  bank,  as 
well  as  initializing  filters,  will  be  discussed. 

2.4.1  Moving  the  Bank 

As  stated  previously,  the  moving-bank  MMAE  is  the  7  subset  of  filters  which 
surround  the  most  likely  correct  filter.  When  parameter  variations  take  place,  the  bank 
will  have  to  move  to  keep  itself  centered  over  the  correct  filter.  This  moving-bank 
concept  is  depicted  in  Figure  1-6.  In  order  to  accomplish  the  bank  move,  three  methods 
(or  logics)  are  presented.  They  are  (1)  Residual  Monitoring,  (2)  Parameter  Position 
Estimate  Monitoring,  and  (3)  Probability  Monitoring. 

2.4. 1.1  Residual  Monitoring 

The  residual  monitoring  method  makes  use  of  the  likelihood  quotient,  Ly(f,), 
which  is  defined  as: 

=  (2.35) 

Note  that  the  likelihood  quotient  has  the  same  quadratic  form  as  the  exponential  term  in 
Equation  (2.27).  Thus  at  each  sample  time,  Lj(ti)  is  computed  for  each  of  the  J  filters. 
Each  of  the  J  likelihood  quotient  values  is  compared  to  a  set  threshold.  Computed  values 
larger  than  the  threshold  are  an  indication  the  bank  needs  to  be  moved  (in  the  direction  of 
the  smaller  values).  If  all  L;  (f,  )  values  are  greater  than  a  set  threshold,  then  the  bank  can 
be  moved  in  the  direction  of  the  smallest  value,  or  the  bank  can  be  expanded  (to  be 
explained  shortly).  The  direction  of  movement  (determined  from  L;(f,)  with  the  smallest 
value)  should  correspond  to  the  parameter  a  j  being  closest  to  the  true  parameter.  Since 
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this  method  relies  on  the  assumption  that  the  most  likely  true  filter  location  is  the  one 
corresponding  to  the  smallest  current  residual  or  in  this  case  the  smallest  current  likelihood 
quotient  value,  then  it  is  inherently  troubled  with  possible  single  large  residual  values 
giving  an  incorrect  indication  of  the  move  direction.  For  this  reason,  the  move  decision 
threshold  must  be  set  high  enough  so  the  bank  will  not  move  in  an  erratic  fashion 
[22:1876]. 

2.4. 1.2  Parameter  Position  Estimate  Monitoring 

This  method  uses  the  parameter  estimate  itself  to  keep  the  moving-bank  centered 
over  the  most  likely  correct  filter  model.  This  is  accomplished  by  comparing  the  current 
parameter  estimate  from  Equation  (2.32)  to  the  bank’s  center.  If  the  difference  is  greater 
than  a  set  threshold,  then  the  bank  can  be  moved  in  the  appropriate  direction  so  as  to  put 
its  center  as  close  as  possible  to  the  current  parameter  estimate. 

In  comparing  this  method  with  the  previous  one,  it  is  seen  that  using  the  parameter 
position  estimate  is  a  more  direct  interpretation  of  the  information  indicating  the  true 
location  of  the  uncertain  parameter.  It  is  also  seen  that  this  method  utilizes  the  entire  past 
measurement  history,  whereas  the  Residual  Monitoring  method  utilizes  only  the  current 
measurement.  Thus,  the  Parameter  Position  Estimate  Monitoring  method  does  not  suffer 
from  possible  erratic  behavior  due  to  single  large  residual  values. 

2.4. 1.3  Probability  Monitoring 

This  method  combines  properties  of  the  previous  two  methods  and  attempts  to 

center  the  bank  over  the  most  likely  filter  model  by  utilizing  the  conditional  probabilities 
determined  from  Equation  (2.26).  The  Pj  (t, )  with  the  highest  value  (greater  than  a  set 

threshold)  is  used  to  determine  the  bank’s  center,  and  the  bank  will  move  in  the  direction 
of  the  filter  with  the  highest  computed  probability.  Thus,  (as  with  Residual  Monitoring) 


2-15 


this  method  uses  residual  information,  except  that  here  the  entire  measurement  history  is 
utilized  as  well.  Then,  (as  with  Position  Monitoring)  this  method  is  not  susceptible  to 
erratic  behavior  due  to  single  large  residual  values. 

2.4.2  Expanding  the  Bank 

When  parameter  variations  are  “small”,  then  the  set  of  filters  can  move  as  a  finely 
discretized  bank.  If,  however,  the  true  parameter  makes  a  large  jump  in  value  in  the 
parameter  space,  the  bank  may  need  to  be  expanded  to  a  coarse  discretization  in  order  to 
encompass  the  true  parameter  within  the  bank.  The  decision  logic  of  whether  or  not  to 
expand  the  bank  is  based  on  the  likelihood  quotient  Equation  (2.35)  (the  other  methods 
rely  on  a(/,)  and  Pj(t()t  both  of  which  are  confined  to  the  nearest  neighbors  to  the 

current  bank  center,  using  fine  discretization).  As  with  Residual  Monitoring,  if  all  L;  (/, ) 

values  exceed  a  set  threshold,  the  bank  is  expanded.  This  too  can  be  prone  to  erratic 
expansions  due  to  single  large  residual  values.  The  threshold  must  be  determined  through 
trial  and  error  for  a  given  system. 

In  addition,  the  bank  can  be  initially  expanded  in  an  effort  to  approximate  the 
location  of  the  true  parameter  within  the  entire  parameter  space.  Once  this  acquisition 
cycle  yields  an  estimate,  a ,  the  bank  can  be  contracted  to  fine  discretization,  centered  on 
a.  Maybeck  and  Hentz  [5:26,  22:92]  stated  that  an  initial  coarse  discretization  may 
provide  a  better  initial  transient  in  parameter  estimation  than  an  initial  fine  discretization. 

2.4.3  Contracting  the  Bank 

After  the  bank  has  been  expanded,  it  will  need  to  be  contracted  once  an  accurate 
estimate  of  the  uncertain  parameter  is  obtained.  The  decision  to  contract  the  bank  is 
generally  controlled  by  a  method  called  the  Parameter  Estimation  Error  Covariance 
Monitoring  technique.  The  covariance  of  the  uncertain  parameter  estimate,  Pa(f,-),  is 
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given  by  Equation  (2.33).  Similar  to  previous  decision  making  methods,  when  some 
scalar  length  variable  associated  with  this  covariance  falls  below  a  set  threshold,  the 
coarsely  discretized  bank  is  contracted  to  a  finely  discretized  bank.  Also,  as  with  previous 
methods,  the  threshold  must  not  be  set  too  “high”  or  too  “low.”  Hentz  [5:69]  found  that 
too  high  a  setting  causes  the  bank  to  contract  prematurely  (without  an  adequate  parameter 
estimate),  and  too  low  a  setting  delays  the  bank  from  contracting  at  the  proper  time. 

Since  Ps  (t, )  is  not  in  general  a  scalar,  in  the  past  there  have  different  methods  for 

utilizing  this  covariance  information.  Hentz  used  either  the  largest  diagonal  element 
[5:64]  or  the  norm  of  P.(ff)  [22:92],  whereas  Kamick  proposed  another  method  involving 

the  probability  of  each  side  of  the  bank  [7:28-29].  The  probability  of  each  side  of  the  bank 
is  given  by 


1  a>  ’  Z‘-»  ^ 
X4 I  a>’  Z|->) 


(2.36) 


With  this  method,  if  the  probability  of  a  side  falls  below  a  given  threshold,  that  side  is 
contracted  “inward”  toward  the  parameter  estimate.  However,  if  probability  of  a  side 
goes  above  another  threshold,  that  side  is  left  alone,  and  the  other  three  sides  are 
contracted  inward.  The  third  case  for  this  method  is,  if  the  sum  of  all  four  sides  are  below 
a  set  threshold,  then  they  are  all  contracted  inward.  Thus,  contracting  the  bank  can  be 
performed  based  on  different  methods  just  as  with  the  bank  move  logic. 

2.4.4  Initialization  of  New  Elemental  filters 

An  initialization  process  must  take  place  when  the  bank  is  moved,  expanded  or 
contracted,  or  first  started  up.  There  are  three  steps  to  this  process.  First,  the  Kalman 
filter  matrices,  Ot,  Bdt,  Gik,  Hk,  Ak,  Dk  and  Kt,  and  the  controller  matrix,  Gck* 
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(for  the  Icth  elemental  filter  and  controller),  must  be  read  in  from  memory  for  the  filters 
that  now  comprise  the  bank.  Each  of  these  matrices  is  constant  for  each  particular  filter 
(based  on  at)  and  thus  can  be  stored  and  need  not  be  computed  on-line.  Second,  each 
filter  which  is  brought  into  use  must  have  its  state  estimate,  xk(ti),  initialized  to  the 
current  state  output  from  the  MMAE,  x(/,).  Third,  the  filter  probability  weights,  pk(t,), 
must  be  initialized.  There  are  different  methods  used  to  initialize  the  probability  weights, 
and  these  are  discussed  below. 

When  the  filters  are  first  initialized,  the  probability  weights  are  equally  divided 
among  the  /  filters  (pj(t0)  =  \/  J  for  all  j).  This  is  justified  since  at  there  is  no  a  priori 

information  about  the  probabilities.  This  will  be  the  same  case  when  the  bank  is  expanded 
or  contracted  with  the  same  justification.  However,  when  the  bank  is  moved,  a  more 
appropriate  probability  weight  distribution  is  warranted  for  the  newly  declared  filters. 
Instead  of  giving  them  all  equal  weighting  (as  in  the  initialization),  the  new  filters  are  given 
an  equal  distribution  of  the  deactivated  probabilities  (the  total  probability  of  the  filters  that 
are  taken  off-line  by  making  the  move).  The  filters  that  are  preserved  during  the  move 
retain  their  previously  computed  probabilities.  Other,  more  complicated,  methods  of 
distributing  the  probabilities  have  been  proposed  and  researched  [17:27-29],  but  the 
method  of  equal  distribution  was  found  to  be  just  as  effective  and  will  therefore  be  used  in 
this  research. 

2.5  Stochastic  Controller  Development 

Now  that  the  estimator  portion  of  the  system  is  developed,  we  will  move  on  to  the 
controller  development.  As  indicated  in  Chapter  1 ,  an  LQG  stochastic  controller  is  to  be 
implemented.  For  the  most  general  case,  the  cost  function  to  be  minimized  is  defined  as 
[21:73]: 
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(2.37) 
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where: 


•  J  —  scalar  cost  to  be  minimized 

•  x(r, )  =  n-dimensional  state  vector 

•  X(r,  )  =  n  x  n-dimensional  state  weighting  matrix 

•  Xf  =  n  x  n-dimensional  final  state  weighting  matrix 

•  u(r,  =  r-dimensional  control  input  vector 

•  U(f, )  =  r  x  r-dimensional  control  weighting  matrix 

•  S(f, )  =  n  x  r-dimensional  cross-weighting  matrix 

•  tN  =  last  time  a  control  is  applied  (held  constant  to  next  sample  period) 

•  tN+l  =  final  time 


Here,  X(f, )  and  Xf  are  real,  symmetric  and  positive  semi-definite;  U(f,)  is  real, 
symmetric  and  positive  definite.  S(/,  )  is  chosen  such  that  [X(f, )  -  S(f,  )U~1(tI  )ST(t, )]  is 

positive  semi-definite  and  such  that  the  augmented  matrix  in  Equation  (2.37)  is  positive 
semi-definite  [21:73].  Values  in  X(f,)  and  Xf  determine  the  relative  penalty  for  not 
controlling  the  states  to  zero,  while  values  in  U(f,  )  determine  the  relative  amount  of 

penalty  for  applying  power  to  the  actuators  to  accomplish  the  controlling  action.  For 
diagonal  X(ff)  and  Xf,  the  larger  the  values  in  these  matrices,  the  more  importance  is 
placed  on  keeping  particular  states  at  zero.  Similarly,  for  diagonal  U(t;),  the  larger  the 
values,  the  more  importance  is  placed  on  keeping  control  deviations  small  on  particular 
actuators.  The  values  that  are  placed  in  these  weighting  matrices  come  from 
experimentation  and  observation  of  the  system  of  interest.  For  some  desired  performance, 
these  matrices  must  be  tuned  in  the  same  fashion  as  Qd,  R  and  P0  in  the  Kalman  filter. 
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This  correspondence  between  estimator  and  controller  is  part  of  a  “duality”  relationship 
and  will  be  discussed  further  shortly  [21:11]. 

With  the  quadratic  cost  defined  above,  the  controller  to  minimize  that  cost  is  given 
by: 

u*(f,)  =  (2.38) 

which  is  obtained  using  certainty  equivalence  with  the  deterministic  case  [21:70]. 

The  controller  gain,  Gc*(f,),  is  determined  by  solving  the  backward  running  Riccati 

recursion  [21:73]: 


G;((i)  =  [U(/i)+BdT«i)Kc((w)Bd((,)]'1[BaT((l)Kt(rl„)<Kiw,»,)+ST(0]  (2.39) 


Kc(l,)=  X(<,)  +  't>T  (I,  )KC  (/,„  )<6(l„, ,  t, ) 

— [BdT<*,  )Kc(<m  )®((m  ,  ( ) + ST  ((,  )]T  G/ft ) 


(2.40) 


solved  backwards  from  the  terminal  condition 


For  a  time-invariant  system  the  control  law  becomes 


(2.41) 


u*(rf)  =  -G/x(0 


(2.42) 


G*  can  now  be  solved  off-line  as: 
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G;  =  [U  +  BaTKcBd]_,[BJTKc<P  +  ST] 


(2.43) 


where  K  is  the  solution  to  the  steady  state  Riccati  equation: 


Kc  =  X  +  <DtKcO  -  [BdTKc<I>  +  STf[U  +  BdTKcBd  ]"'[BdTKe<I>  +  ST]  (2.44) 


The  time-invariant  condition  is  welcomed  for  high  dimensional  systems  (for  decreased 
computational  load).  This  condition  is  physically  motivated  when  the  initial  Kalman  filter 
transients  and  the  final  controller  gain  transients  are  both  negligible.  This  will  be  the  case 
if  the  two  transient  time  intervals  are  small  compared  to  the  time  interval  that  control  is 
applied. 

A  note  is  made  here  concerning  the  “duality”  relationship  between  the  Kalman 
filter  Riccati  recursion  Equations  (2.18),  (2.19)  and  (2.21),  and  the  LQG  controller  Riccati 
recursion  Equations  (2.39)  and  (2.40).  This  duality  relationship  will  be  shown  explicitly 
as  follows.  First,  we  substitute  Equations  (2.19)  and  (2.21)  into  Equation  (2.18)  (note  the 
short  notation  used  for  compactness): 


p-  =  o p-j<&T  -  oi^:1ht[h^:1ht + r]"1  h^:,ot  +  GdQdGdT  (2.45) 


Then  we  substitute  Equation  (2.39)  into  Equation  (2.40). 


Kc;  =  X  +  ®TKCM®-[BaTKcw®+ST]T 
•  [U  +  BdTKCMBa]>aTKtw®  +  ST] 


(2.46) 
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Now,  by  comparing  Equations  (2.45)  and  (2.46),  the  following  variable  relationships  can 
be  made: 


X<-»GdQdGj 

<I>  4>t 

Bd<-*HT 
U<->R 

With  the  corresponding  relationships  just  noted,  it  is  seen  that  the  Kalman  filter  Riccati 
recursion  is  of  the  exact  same  form  as  the  LQG  Riccati  recursion.  This  duality 
relationship  is  worth  noting  because,  as  a  result  of  the  two  Riccati  equations  being  of  the 
same  form,  they  can  both  be  solved  for  their  steady  state  solutions  by  the  same  routine. 

One  last  simplification  is  made  concerning  the  cross  weighting  matrix,  S,  in  the 
quadratic  cost  function.  Equation  (2.37).  This  term  has  been  shown  in  previous  research 
[1 1]  to  have  negligible  impact  on  performance  and  is  thus  neglected. 

2.6  MMAC 

In  parallel  with  the  Kalman  filter/MMAE  development  is  the  LQG  control/Multiple 
Model  Adaptive  Control  (MMAC)  development.  Now  that  the  control  equations  have 
been  presented,  the  methods  by  which  the  control  is  implemented  are  presented.  Just  as 
there  are  several  methods  of  implementing  the  state  and  parameter  estimation  algorithms, 
there  are  various  methods  of  implementing  the  controller.  The  following  six  methods  are 
presented  in  this  section:  (1)  MMAC  Control,  (2)  “Modified”  MMAC  Control,  (3) 
MMAC  with  MAP  Control,  (4)  Single  Fixed-Gain  Control,  (5)  Single  Changeable-Gain 
Control,  and  (6)  “Modified”  Single  Changeable-Gain  Control. 
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2.6.1  MMAC  Control 

The  basic  structure  of  the  MMAC  is  shown  in  Figure  1-2.  This  is  the  basis  from 
which  the  other  methods  are  derived.  With  the  basic  MMAC  algorithm,  an  elemental 
controller  is  cascaded  with  each  elemental  Kalman  filter.  This  structure  parallels  that  of 
the  MMAE  in  that  the  output  of  each  elemental  controller  is  probabilistically  weighted  and 
summed  to  provide  the  optimal  control,  i.e., 

K 

»'«,)= <2-47) 
k= 1 

where, 

u/(r,.)  =  -G/(at)xt(0  (2.48) 

Note  that  each  elemental  controller  gain  is  based  on  the  corresponding  assumed  parameter 
value,  a*,  for  that  filter  model.  This  form  of  Multiple  Model  Adaptive  Control  is  a  natural 

extension  of  the  Bayesian  Multiple  Model  Adaptive  Estimator. 

2.6.2  “Modified”  MMAC  Control 

The  “modification”  in  this  method  is  an  extension  to  the  ad  hoc  technique  of  lower 
bounding  the  probabilities,  pk(t{),  presented  in  Section  2.3.2.  With  this  method,  after 

each  of  the  K  control  vectors  are  formed,  the  associated  probabilities  are  compared  to  a 
set  threshold,  at  a  numerical  value  higher  than  the  lower  bound  of  that  previous  section.  If 
any  pk  (ti )  is  found  to  be  below  this  threshold,  the  corresponding  kth  control  vector, 
),  is  discounted  and  not  used  to  form  the  probabilistically  weighted  control  vector. 
The  probabilities  are  rescaled  so  their  sum  is  one,  and  the  remaining  u*  (/, )  are  used  to 
form  the  control  vector  per  Equation  (2.47). 
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The  use  of  the  lower  bounds  described  in  Section  2.3.2  is  still  important  to  prevent 
the  “lock-out”  phenomenon  associated  with  a  particular  pk  {ti )  value  being  computed  as 
essentially  zero.  However,  if  that  pt(f;)  is  artificially  lower  bounded,  the  control 

computed  via  Equation  (2.47)  will  include  a  nonzero  weighting  on  the  corresponding 
inappropriate  control,  u k(tt).  In  this  particular  problem,  it  is  readily  possible  that  putting 
any  nonzero  weight  on  that  u k(tt)  could  drastically  degrade  performance  and  could  even 

cause  closed-loop  instability.  Thus,  the  lower  bounding  method  is  incorporated  to 
maintain  viable  pk  (ti )  computations,  but  any  (r, )  associated  with  a  pk  (f, )  that  has  been 

artificially  bounded  is  deemed  to  be  unsuitable  for  use  in  control  generation.  Using  a 
threshold  somewhat  higher  than  the  lower  bound  value  allows  other  Ut(ff)’s  that  are  poor 
candidate  controls  (as  indicated  by  their  low  pk(tt)  values)  also  to  be  eliminated  from  the 
control  calculation. 

2.6.3  MMAC  with  MAP  Control 

The  method  corresponds  directly  with  the  MAP  method  of  state  estimation 
(Section  2.3.2);  the  control  is  based  solely  on  the  elemental  LQG  controller  with  the 
highest  associated  probability. 

2.6.4  Single  Fixed-Gain  Control 

This  method  utilizes  the  fact  that  full-state  feedback  controllers  are  inherently 
robust  and  determines  the  control  input  based  on  a  single  nominal  parameter  value  and  a 
single  state  estimate  as  output  from  the  MMAE  [5:39-40,  22:93].  Figure  2-1  illustrates 
this  concept,  and  the  control  law  is  given  as: 

u«,.)  =  -G/(a„Jx(0  (2.49) 


Hentz  noted  that  determining  a  good  value  for  may  not  be  a  trivial  task  [5:40]. 
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Figure  2- 1 .  Single  Fixed-Gain  Controller 


2.6.5  Single  Changeable-Gain  Control 

This  method  is  similar  to  the  previous  in  that  it,  too,  uses  the  single  state  estimate 
as  output  from  the  MMAE.  However,  the  controller  gain  is  computed  based  on  the 
parameter  estimate  as  output  from  the  MMAE  instead  of  a  single  nominal  parameter 
value.  This  idea  is  illustrated  in  Figure  2-2,  and  the  control  law  is  given  by: 


u*(f, )  =  -Gc*[a(f,)]x(f,+ )  (2.50) 

This  requires  that  gain  values  be  stored  for  each  of  the  K  parameters.  Then,  Gc*[a(^.)]  is 
computed  by  interpolation  from  the  stored  values[5:37-38]. 
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2.6.6  “ Modified ”  Single  Changeable-Gain  Control 

This  method  is  similar  to  the  previous  method  in  that  it  also  uses  the  parameter 
estimate  from  the  MMAE,  but  the  state  estimate  comes  from  a  separate  Kalman  filter  (not 
part  of  the  moving  bank  of  filters).  This  controller  structure  is  illustrated  in  Figure  2-3, 
and  the  control  law  is  defined  the  same  as  in  the  previous  method  with  x(/f+)  from  the 
MMAE  replaced  by  x(t*)  from  the  separate  Kalman  filter.  This  form  of  control  can  be 
useful  if  a (t/)  is  between  two  discrete  parameter  values  for  which  elemental  filters  are 
designed.  A  single  filter  based  on  a(t( )  might  provide  substantially  better  performance 
than  either  of  those  two  elemental  filters  (or  any  linear  combination  of  those  filters);  this 
form  of  algorithm  allows  a  coarser  discretization  of  the  parameter  space  [18].  This 
method  of  control  also  reduces  the  possibility  of  control  inputs  based  on  an 
underestimated  undamped  natural  frequency,  which  previous  research  has  shown  to  cause 
instabilities  in  this  type  of  system  [5,  22:95,  30]. 
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Figure  2-3.  Modified  Single  Changeable-Gain  Controller 


2.7  Dithering  Techniques 

As  indicated  in  Chapter  1,  this  research  will  apply  dithering  techniques  in  order  to 
enhance  parameter  identification.  Within  the  area  of  dithering,  there  are  two  basic 
methods:  purposeful  excitation  designed  for  enhancement  of  parameter  identification  and 
purposeful  commands  to  accomplish  slewing  of  the  space  structure.  Both  types  of 
commands  are  referred  to  as  “open  loop”  control,  since  each  are  precomputable  functions 
of  time  and  not  of  states  or  measurements.  If  an  appropriate  excitation  frequency  is  used, 
the  system  should  be  excited  enough  to  enhance  parameter  identification.  Likewise, 
giving  the  structure  a  slewing  motion  will  cause  a  certain  amount  of  natural  system 
excitation  to  occur,  which  may  also  enhance  parameter  identification. 

This  section  will  discuss  the  different  types  of  purposeful  excitation  dithering  to  be 
used  in  this  research.  The  rigid-body  states  and  how  they  are  incorporated  into  the 
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structural  and  control  matrices,  F  and  B,  for  use  with  purposeful  (slew)  commands,  will 
be  presented  in  Chapter  3. 

2. 7. 1  Sinusoidal  Dither 

Three  sinusoidal  dither  wave  forms  will  be  considered  in  this  research.  They  are: 
(1)  sine  wave,  (2)  swept  frequency  sine  wave,  and  (3)  square  wave.  The  sine  wave  would 
be  a  logical  choice  if  the  system’s  identification-enhancing  frequency  were  known  a 
priori.  This  optimum  frequency  is  one  which  (when  used  to  excite  the  system)  should 
allow  the  best  parameter  identification  to  be  made.  If,  however,  this  optimum  frequency  is 
not  known,  a  method  of  using  a  sine  wave  while  slowly  changing  the  frequency  across  a 
particular  range  (swept  frequency  sine  wave)  could  be  employed  [18].  This  method  would 
allow  the  system’s  identification-enhancing  frequency  to  “show  itself’,  and  the  dither 
could  then  be  switched  to  a  sine  wave  at  that  particular  frequency  (or  be  swept  in  some 
small  frequency  band  about  that  frequency). 

An  alternative  to  using  a  pure  sine  wave  dither  is  to  use  a  square  wave.  With  a 
square  wave,  not  only  is  the  fundamental  frequency  of  the  wave  induced  into  the  system, 
but  additional  harmonics  (high  frequency  content)  are  also  induced  [18].  These  additional 
harmonics  may  help  parameter  identification  when  the  optimal  dither  frequency  is  not 
known. 

When  implementing  sinusoidal  type  signals  in  a  digital  computer  simulation,  certain 
limitations  exist  and  are  discussed  in  the  next  section. 

2. 7.2  Dither  Frequency 

In  a  discrete-time  system  which  is  propagated  forward  at  fixed  intervals, 
/,  —  =A t,  Shannon’s  sampling  theorem  will  apply  when  computing  an  input  dither 
signal  in  the  form  of  a  sin(cor(fl ))  function.  In  general,  Shannon’s  sampling  theorem 
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requires  that  the  sampling  rate  be  at  least  twice  the  highest  frequency  of  interest  in  the 
signal  being  sampled.  This  limit  is  also  known  as  the  Nyquist  rate  [27:  81].  For  our 
application  of  computing  a  dither  signal  in  a  simulation  which  contains  a  particular  given 
sampling  rate,  the  limitation  is  placed  on  the  frequency  of  the  dither  signal  which  can  be 
produced.  Thus,  the  maximum  frequency  in  a  sinusoid  function  that  can  be  correctly 
produced  in  a  fixed  sampling  rate  simulation  is: 


(2.51) 


where  /  is  highest  frequency  or  Nyquist  rate  (in  Hertz)  with  which  the  sin(©r(r(  )) 
computation  will  produce  a  useful  result.  (Note  that  (0  =  2nf,  if  (0  is  in  rad/sec.)  Thus, 
computing  sin((or(/1 ))  with  the  limitation  set  in  Equation  (2.51)  is  analogous  to  sampling 
some  continuous  signal  of  frequency,  /,  with  the  Nyquist  sampling  rate  requirement 
fs  >  2/,  where  fs  is  the  sampling  rate.  For  the  system  used  in  this  research 
At  =  O.Olsec,  and  therefore,  /  <  50 Hz  or  CO  <  314.16rad/  sec.  The  importance  of  this 
maximum  frequency  will  be  seen  in  Chapter  4  when  the  actual  dither  frequencies  are 
chosen. 

One  additional  note  concerns  the  actual  application  and  computation  of  the  sine 
wave  and  square  wave  dither  inputs.  Dither  is  applied  via  the  control  variable,  u(ff ),  by: 


«(/,)  =  Asin(co/()  (sine  wave)  (2.52) 

u(tt )  =  A  sgn(sin(co/. ))  (square  wave)  (2.53) 


where  A  is  the  amplitude  of  the  wave,  and  the  signum  function  “sgn”  merely  returns  a  +1 
or  a  -1,  depending  on  the  polarity  of  the  sine  wave. 
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2.7.3  Wide  Band  Noise  Dither 


An  alternative  to  using  a  sinusoid  type  of  dither  signal  is  to  use  wide  band  noise. 
With  wide  band  noise,  a  shaping  filter  is  designed  for  the  purpose  of  providing  a  signal 
that  covers  a  specified  range  of  frequencies  at  a  specified  strength  or  power  spectral 
density  (PSD)  value  (over  all  inclusive  frequencies).  With  such  a  shaping  filter,  the  dither 
input  to  the  system  can  provide  a  very  rich  frequency  content,  since  a  broad  range  of 
frequencies  is  used. 

For  this  research  a  simple  second-order  shaping  filter  with  two  real  poles  will  be 
used.  Figure  2-4  shows  the  transfer  function  and  equivalent  block  diagram  design  for  this 
shaping  filter.  The  PSD  plot  which  represents  the  output  of  this  transfer  function  and 
block  diagram  is  shown  as  Figure  2-5.  The  corresponding  state-space  equations  are  as 
follows: 


*j(0 

'-A 

0  Txj(f)' 

4» 

w  (O' 

x2(/> 

l 

~fi_  .MO. 

i 

_  0 

n(/)  =  K*,(»)  =  k(x,(<)-/2xj»)  (2.55) 


where  W  is  zero  mean,  unity  strength  white  noise,  and  K  is  an  adjustable  gain.  fx  and  /2 
are  the  lower  and  upper  frequency  ranges  of  interest.  Equation  (2.54)  can  be  written  as 
the  continuous-time  matrix  differential  equation: 

X(f)  =  Fx(f)  +  W(0  (2.56) 
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Transfer  Function  Form 


Figure  2-4.  Block  Diagram  for  Wide  Band  Noise 


Figure  2-5.  PSD  Plot  for  Wide  Band  Noise 
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For  discrete-time  implementation,  Equation  (2.S6)  is  written  as: 


X(t, )  =  <!>(/,  -  )X(f,_1 )  +  Wd  (fM )  (2.57) 

where  0(f,  the  state  transition  matrix  associated  with  F  and  the  statistics  for  Wd 

can  be  computed  according  to  Equations  (2.1 1),  (2.12)  and  (2.13).  For  constant  sampling 
period  {ti  -  f(_,  =  const  =  At),  the  state  transition  matrix  is  simply: 


<S>(At)  = 


-Mi 


0 


Si-fi 


(2.58) 


The  entire  gain  of  the  output  of  the  shaping  filter  is  controlled  by  K  in  Equation  (2.55) 
and  will  be  adjusted  by  trial  and  error. 

2.8  Mathematical  Modeling  Methods 

This  section  will  present  the  two  coordinate  forms,  physical  and  modal.  The 
transformation  from  the  physical  form  to  the  more  desirable  modal  form  will  also  be 
presented,  along  with  a  technique  for  producing  a  reduced-order  model  based  on  the 
modal  coordinate  form. 

Many  systems  (including  the  one  with  this  research)  contain  a  high-dimension  truth 
model.  This  is  very  useful  in  computer  simulations  in  order  to  represent  the  “real”  world 
accurately.  However,  many  on-line  applications  cannot  perform  with  such  high  state 
dimensions.  Thus,  a  reduced-order  model  as  the  basis  of  filter  and  controller  designs  is 
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necessary  which  will  satisfy  the  on-line  computer  limitations  and  which  will  also  yield 
satisfactory  performance. 

2. 7. 1  Physical  Coordinate  Form 

The  governing  second-order  differential  equation  describing  the  flexible  structure 
vibrations  is  given  by  [7:39, 13:4,  32]: 

Mr(/)  +  Cr(0+Kr(0  =  F,(u,/)+F2(0  (2.59) 

where: 

•  r(f)  =  n-dimensional  vector  representing  the  structure's  physical  position 

•  M  =  n-by-n  constant  mass  matrix 

•  C  =  n-by-n  constant  damping  matrix 

•  K  =  n-by-n  constant  stiffness  matrix 

•  F,(t)  =  r-dimensional  deterministic  control  inputs 

•  F2(r)  =  r-dimensional  disturbances  and  unmodeled  control  inputs 

With  the  assumptions  that  the  deterministic  control  and  noises  enter  the  system  linearly, 
and  that  the  external  disturbances,  F2(r),  can  be  modeled  as  white  Gaussian  noises,  then 

Equation  (2.59)  becomes  [7:40,  31:4]: 

Mr  (0 + cr(0 + Kr(0  =  -bu(t)  -  gw(r)  (2.60) 

where: 

•  u(r)  =  r-dimensional  vector  of  actuator  inputs 

•  b  =  n-by-r  control  input  matrix  identifying  position  of  actuators  and  relationship 
between  actuators  and  controlled  variables 

•  W(r)  =  5-dimensional  vector  representing  the  dynamics  driving  noise 


2-33 


•  g  =  n-by-s  noise  input  matrix  identifying  position  of  input  disturbances  and 
relationship  between  the  dynamics  driving  noise  and  the  controlled  variables 

Equation  (2.60)  can  be  written  as  the  stochastic  differential  equation  in  state-space  form 
[7:40,  13:4]: 


X(/)  =  Fx(f)  +  Bu(f)  +  Gw(f)  (2.61) 


Note  that  this  equation  is  of  the  same  form  as  the  system  described  by  Equation  (2.1)  in 
Section  2.2  with  the  exception  of  constant  matrices.  The  common  state  vector 
representation  is  given  as: 


x(0  = 


r(0 


r(t) 


2nxl 


(2.62) 


and  the  form  of  the  constant  system  matrices  is  [2:3-20, 7:41]: 


F  = 


-M-'C...  -M-'K...1, 


*  “mb 

0 

nxn  J2nx2n 


(2.63) 


B  = 


nxr  J2 nxr 


(2.64) 


G  - 


0. 


nxs 


12/1X5 


(2.65) 


/or  the  state  vector  described  in  Equation  (2.62),  the  discrete-time  measurement 
model  of  position  and  velocity  is  given  by  [2:3-20,  7:41]: 


where: 


0 

H, 


x(/,) 


P  Jmx2n 


+v(0 


(2.66) 


•  m  =  number  of  measurements 

•  v(/,)  =  m-dimensional  measurement  uncertainty  modeled  as  a  discrete-time 

white  Gaussian  noise  of  covariance  R(f,) 

•  Hp  =  (m/2)-by-n  position  measurement  matrix  in  physical  coordinates 

•  Hv  =  (m/2)-by-n  velocity  measurement  matrix  in  physical  coordinates 


Note  that,  for  simplicity,  an  equal  number  of  position  and  velocity  measurements  is 
assumed.  The  actual  measurement  matrix  may  vary,  depending  on  available 
measurements,  and  does  not  affect  the  structure  of  the  system  matrices  of  Equations 
(2.63)  through  (2.65). 

From  this  development  it  is  seen  that  the  physical  coordinate  form  produces  highly 
coupled  system  equations.  With  this  coupled  format,  important  system  characteristics  may 
not  be  easily  identified,  and  thus,  a  more  desirable  form  is  needed.  The  next  section  will 
discuss  the  transformation  from  physical  to  modal  coordinate  form. 

2.7.2  Modal  Coordinate  Form 

Transforming  the  equations  to  modal  form  decouples  the  modes  and  makes  the 
identification  process  simpler.  For  the  purposes  of  this  research,  the  damping  matrix,  C, 
is  assumed  to  be  a  linear  combination  of  the  mass  and  stiffness  matrices  [13:4]: 
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C  =  cxM  +  pK 


(2.67) 


However,  the  calculation  of  (X  and  (3  is  not  necessary  when  transforming  to  the  modal 
coordinate  form,  as  will  be  seen  shortly.  Given  the  new  modal  coordinate  vector,  f ,  the 
transformation  from  physical  form  to  modal  form  is  described  by  [13:5]: 


r  =  TV  (2.68) 

where  T  is  a  n-by-n  transformation  matrix  determined  from  the  system  eigenvectors  by 
solving  [13:5,  32]: 

C02Mr  =  K  T  (2.69) 

where  the  values  for  (0  that  satisfy  this  equation  are  referred  to  as  the  natural  or  modal 
frequencies. 

Now,  using  the  transformation  Equation  (2.68)  to  operate  on  the  original  ‘■ystem 
Equation  (2.61),  the  resulting  state  space  equation  is  given  by  [13:5]: 

x(t)  =  Fx(0  +  Bu(r)  +  Gw(f )  (2.70) 


where  the  transformed  state  vector  from  Equation  (2.62)  is  now  defined  as  [13:5]: 


X(t)  = 


(2.71) 
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and  the  transformed  matrices  from  Equation  (2.61)  as  applied  to  Equation  (2.60)  are 
defined  as  [2:3-22,  13:5]: 


F  = 


i  o 


2nx2n 


B  = 


r'M'b 
0 


J2  nxr 


G  = 


Tx  M^g 
0 


J2  nxs 


(2.72) 


(2.73) 


(2.74) 


Since  the  modal  form  provides  for  independent  equations,  the  resulting  modal 
vectors,  i.e.,  eigenvectors  of  F,  are  orthogonal.  This  fact  plus  the  following  relationships 
[7:44,  32:1769]: 


=  [-20,-]  (2.75) 

-T-'M-'KT  =  [-©?]  (2.76) 


allows  the  dynamics  matrix  to  be  written  as  follows: 


[-2 Oi]  H>,2J 


I 


0 


\2nx2n 


(2.77) 


where  the  matrix  is  now  block  diagonal  in  terms  of  the  undamped  natural  frequency  and 
the  damping  ratio  of  the  i-th  mode.  In  Equations  (2.75)  through  (2.77),  the  notation 
“[a,]”  denotes  a  diagonal  matrix  with  a /  through  a„  running  down  the  diagonal.  Note  that. 
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because  of  the  relationship  in  Equation  (2.7S),  the  damping  matrix  C  does  not  need  to  be 
determined  and  thus  the  a  and  j3  in  Equation  (2.67)  need  not  be  known.  The  transformed 
measurement  equation  from  Equation  (2.66)  may  be  written  as  [2:3-22]: 


Z(/,)  = 


THvr 

0 


0 

H  PT 


\mx2n 


+  v(f,) 


(2.78) 


From  this  development  it  can  be  seen  that  the  modal  coordinate  form  does  provide 
a  decoupled  set  of  equations  from  which  to  work.  This  decoupled  form  (modal  form) 
allows  direct  insight  into  the  physical  system  parameters  (natural  frequency  and  damping 
ratio),  unlike  with  the  physical  coordinate  form. 

The  next  section  presents  the  development  of  the  Modal  Reduction  Technique  as 
presented  by  Kokotovic  [9]  and  Gustafson  [2].  There  are  many  other  methods  of 
providing  a  reduced-order  model  from  a  truth  model,  but  Gustafson  obtained  the  best 
results  using  this  method  [2:5-30],  and  so  this  technique  will  be  pursued  in  this  research. 

2.7.3  Modal  Reduction  Technique 

This  state  reduction  technique  will  eliminate  the  higher  frequency  bending  modes 
of  the  system.  The  concept  behind  forming  a  reduced-order  model  by  eliminating  these 
higher  frequency  modes  (states)  comes  from  the  assumption  that  the  higher  frequency 
modes  will  become  quiescent  in  a  negligibly  small  amount  of  time.  This  assumption  comes 
from  the  physical  shape  of  the  SPICE  structure,  which  is  two  bulkheads  connected  by 
tripod  legs.  With  this  type  of  structure,  it  is  suspected  that  low  frequency  mode  effects 
will  dominate. 
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The  state-space  system  model  from  Equation  (2.61)  can  be  partitioned  as  follows 
[9:123]: 


where  the  system  is  driven  by  deterministic  controls,  u(t),  and  zero-mean,  white  Gaussian 
noise  W(r)  of  strength  Q(t).  The  upper  partition,  X, (/),  corresponds  to  the  low 
frequency  modes  to  be  maintained,  and  the  lower  partition,  X2(r),  corresponds  to  the  high 
frequency  modes  to  be  removed. 

Assuming  instantaneous  steady  state  in  the  higher  modes  (X2(r)  =  0),  these 
modes  can  be  eliminated  with  negligible  impact  to  overall  system  performance.  The  lower 
partition  differential  equation  is  set  equal  to  zero  by  [9: 123]: 

X2  (0  =  F21X,  (0  +  F22X2  (0  +  B2u(r)  +  G2W(0  =  0  (2.80) 

Since  F21  and  F22  are  square  matrices  and  assuming  F2"2  exists,  X2(f)  can  now  be  written 
in  terms  of  X,  ( t )  and  system  inputs: 

X2(0  =  -F22  [F21x,  (r) + B2u(r) + G2  w(0]  (2.81) 

Substituting  Equations  (2.80)  and  (2.81)  into  Equation  (2.79)  results  in  [9:123]: 

^i(f)  =  [Fh  —  F12F22 F2i]xj(r)  +  [Bj  —  F12F22B2]u(r)  +  ^Gi  —  F12F22G2]w(t)  (2.82) 
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Analogous  to  other  developments  in  this  Chapter,  the  equivalent  discrete-time 
model  will  be  developed.  The  equivalent  discrete-time  model  of  Equation  (2.79)  is  given 
as  [2:3-24]: 


X(',+i)  = 


*n 

'x,(/f)" 

Bdl 

U  (*,)  + 

Gdl 

*21 

*22. 

_x2(0. 

T 

Bd2. 

_Gd2_ 

Wd(0 


(2.83) 


Similar  to  the  steady-state  assumption  for  the  continuous-time  case  (X2(f)  =  0),  the 
discrete-time  steady-state  assumption  becomes  (x2(f1+1)  =  X2  (/,)).  Thus,  the  higher- 
order  modes  of  Equation  (2.83)  are  represented  as  [2:3-24]: 


02,X,  (t, )  +  [<*>22  -  I]x2  (/, )  +  Bd  2u(r, )  +  Gd2  Wd  (t, )  =  0  (2.84) 

x2(r,)  =  -[0>22  -  lY^^M+B^uit^+G^yN^)]  (2.85) 

In  order  to  evaluate  these  discrete-time  matrices  from  the  continuous-time  matrices,  a 
first-order  approximation  method  is  used  as  follows  [2:3-24,  19:172,  357]: 


0>22  =I  +  F22Ar  (2.86) 

*21  =  F21Af  (2.87) 

Bd2  =  B2Ar  (2.88) 

Gd2=G2Ar  (2.89) 

Qd=Q/A t  (2.90) 


This  last  equation  comes  from  the  first-order  approximation  of 
Gd2QdGj2  =  G2Q2G2  A t  and  Gd2  =  G2Ar  [2:3-25, 18]. 
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Substituting  Equations  (2.86)  through  (2.90)  into  Equation  (2.85)  results  in: 


X2(0  =  -[Fjj  (/, )  +  B2  Am(f, )  +  G2  A/Wd  (/, )]  (2.91) 


where  W d(0  has  covariance  Qd.  Eliminating  A t  from  this  expression  results  in: 

x2  (*, )  =  -F2‘2  [F2iXi(/(  )  +  B2u(/,.  )  +  G2Wd  (r,. )]  (2.92) 

This  is  the  same  result  obtained  in  the  continuous-time  case,  Equation  (2.81),  with  the 
exception  of  Wd(0  replacing  W(/)-  Then,  substituting  Equation  (2.92)  into  the  discrete- 

time  measurement  equation: 

Z(',)  =  [H,  H2] 

and  expanding,  yields  [2:3-25]: 

2(0  =  [Hj  - H2F^1F21]x1(f|.)-H2F2^,[B2u(f|.)+G2Wd(fI.)]  +  V(0  (2.94) 

The  second  term  in  this  measurement  equation,  -  H2F2~2!  [B2u(f, )  +  G2Wd  (i, )] ,  is  a  direct 

feedthrough  term  created  by  the  order  reduction  [9: 123]. 

One  item  yet  to  be  determined  in  this  modal  reduction  method  is  the  break  between 
the  low  and  high  frequency  modes.  This  is  accomplished  by  examining  an  ordered  list  of 
the  modal  frequencies  and  deciding  what  constitutes  a  natural  break  or  jump  in  modal 
frequency.  The  actual  break-point  made  comes  from  physical  insight  and  desired 
computation  load. 


*i(0 

*2  (01 


+  v(0 


(2.93) 
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With  the  break-point  known,  the  system  dynamics  matrix  from  Equation  (2.77)  is 
separated  into  an  augmented  form  as  [9: 123]: 


[-2$,©,]  [- 

-©?] 

0 

0 

I 

0 

0 

0 

0 

0 

[-2C2co2]  [- 

-<*>2] 

0 

0 

I 

0 

(2.95) 


The  low  frequency  modes  to  be  maintained  are  represented  by  the  upper  left  quadrant,  and 
the  higher  frequency  modes  to  be  removed  are  represented  by  the  lower  right  quadrant. 
These  two  quadrants  correspond  to  the  Fn  and  F22  partitions  in  Equation  (2.79)  with  the 
off-diagonal  blocks,  F12  and  f21,  set  equal  to  zero.  Substituting  this  information  into 
Equations  (2.82)  and  (2.94)  yields  [9:123]: 

i,«)  =  F|li1(»)+Blu«)+GlW,(») 

(2*0) 

=  Frx,(0 + Bru(r)  +  Grwr(f) 

z :(*,)  =  HjijC/,)-  H^fB^) +  G2wd(f,.)l  +  vr(r,.) 

L  J  (2.97) 

=  Hrx,(f, )  -  Duu(r)  +  Dwwd(r,)  +  Vr  (ti ) 

where  the  effect  of  F12  =  0  and  F21  =  0  is  seen  explicitly,  and  where  the  subscript,  r, 
denotes  “reduced-order.”  Thus,  the  direct  feedthrough  matrices,  Du  and  Dw ,  allow  the 
effects  of  control  and  dynamics  driving  noise  inputs  to  be  “seen”  in  the  reduced-order 
measurement. 
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2.9  Summary 

In  this  chapter  the  fundamental  concepts  concerning  the  Kalman  filter,  the  state 
and  parameter  estimator  (the  MMAE),  and  the  stochastic  controller  (the  MMAC)  were 
presented.  For  the  MMAE  and  MMAC,  various  algorithms  and  enhancements  were 
presented  that  have  been  studied  in  the  past.  Shown  below  is  a  summary  of  the 
assumptions  and  particular  estimator/controller  options  (particular  algorithms)  which  are 
to  be  used  in  this  research  effort. 

-  Kalman  filter:  Based  on  linear,  time-invariant  system  model  driven  by 

stationary  noise. 

-  LQG  controller:  Constant  state  and  control  weighting  matrices,  no  cross¬ 

weighting  matrix. 

-  MMAE:  ME/A  and  ME/I  will  be  considered. 

Move  logic  will  be  Parameter  Position  Estimate  Monitoring 

-  MMAC:  “Modified”  MMAC  will  be  used. 

Additionally,  presented  in  this  chapter  were  the  dithering  techniques  to  be  used  for 
enhancing  parameter  identification  and  the  mathematical  modeling  for  developing  reduced- 
order  filter  models.  The  dither  inputs  are  to  consist  of  (1)  square  wave,  (2)  sine  wave,  (3) 
swept  sine  wave,  and  (4)  wide  band  noise.  The  mathematical  modeling  consists  of  (1)  the 
Physical  Coordinate  Form,  (2)  the  Modal  Coordinate  Form,  and  (3)  the  Modal  Reduction 
Technique.  Chapter  3  will  present  the  actual  system  structural  model  development  to  be 
used  in  this  research. 
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3.  System  Development 


3.1  Introduction 

This  chapter  presents  a  summary  of  the  SPICE  2  system  model  originally 
developed  by  Gustafson  [2:  Chapter  3].  The  physical  structure  of  the  SPICE  2  model  is 
described,  followed  by  a  description  of  the  mathematical  models  used  to  simulate  the 
actual  structure  and  its  disturbance  inputs.  The  truth  and  filter  models  used  in  this 
research  are  also  discussed. 

The  addition  of  rigid-body  states  is  presented  here  as  well.  These  are  used  to 
simulate  planar  motion  of  the  structure  so  as  to  induce  naturally  occurring  structural 
vibrations.  This  is  an  addition  to  Gustafson’s  development  and  is  used  to  complement  the 
research  of  applying  dithering  techniques  for  enhancing  uncertain  parameter  identification. 


3.2  SPICE  Structure 

The  structure  model  used  in  this  research  is  model  2  of  the  SPace  Integrated 
Controls  Experiment  (SPICE)  [14:  Chapter  3],  The  actual  structure  resides  at  Phillips 
Laboratory,  Kirtland  AFB,  New  Mexico  and  has  subsequently  been  revised  up  to  model  5. 
For  the  purposes  of  this  research,  however,  the  SPICE  2  model  will  suffice,  since  this 
research  effort  is  devoted  to  applying  and  proving  concepts  within  the  area  of  adaptive 
filters  and  controllers. 
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3.2.1  Physical  Description 

The  SPICE  structure  is  composed  of  three  major  sections,  as  illustrated  in  Figure 
3-1.  The  hexagonal  base  or  bulkhead  is  6.19  meters  in  diameter.  This  base  supports  the 
entire  structure  and  also  contains  the  primary  mirror  assembly  mounted  on  its  top.  Three 
tripod  legs  connect  the  bulkhead  to  the  top  portion,  which  houses  the  Secondary  Mirror 
Assembly  (SMA),  which  is  1.32  meters  in  diameter.  The  overall  height  of  the  structure  is 
8.44  meters.  The  X-Y  axes  are  centered  and  in  the  plane  of  the  structure’s  base,  with  the 
Y-axis  pointing  towards  tripod  leg  number  1.  The  Z-axis  points  straight  up  the  middle  of 
the  structure  and  is  the  Line  of  Sight  (LOS)  axis.  The  12  Proof  Mass  Actuators  (PMAs) 
located  on  the  tripod  legs  contain  their  own  local  coordinate  frames,  while  the  6  PMAs 
located  at  the  base  are  aligned  along  the  Z-axis. 

Keeping  the  secondary  mirror  assembly  aligned  with  the  primary  mirror  assembly 
is  the  main  concern  of  the  SPICE  research.  An  exaggerated  non-alignment  is  shown  in 
Figure  3-2.  For  this  research,  pure  torsional  bending  (about  the  Z-axis)  is  not  considered, 
since  it  has  only  minor  impact  on  the  line-of-sight  tracking  error. 


3.2.2  Actuators  and  Sensors 

The  actuators  are  the  Proof  Mass  Actuator  (PMA)  type  and  are  used  to  quell  any 
vibration  of  the  tripod  legs  or  the  bulkhead.  The  sensors  used  are  accelerometers  and  are 
collocated  with  the  PMAs.  The  accelerometers  are  used  to  measure  the  relative  bending 
motion  of  the  structure.  Twelve  PMAs  and  accelerometers  are  located  on  the  tripod  legs, 
four  on  each  leg.  On  each  tripod  leg,  two  PMAs  and  accelerometers  are  located  at  both 
the  one-third  and  two-third  length  positions.  Six  PMAs  and  accelerometers  are  located  at 
the  base  of  the  structure.  Thus,  the  structure  contains  18  PMAs  and  18  accelerometers. 
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Figure  3-1.  SPICE  Structure  [14:  HI-14] 


Figure  3-2.  Flexible  SPICE  Structure 


3.3  System  Mathematical  Model  Description 

The  mathematical  model  of  each  component  of  the  SPICE  2  system  is  presented 
here.  The  overall  block  diagram  of  the  system  model  is  shown  in  Figure  3-3.  The 
Disturbance,  PMA  and  Accelerometer  blocks  all  contain  shaping  filters  and  state-space 
models  for  their  respective  functions.  The  LOS  block  is  merely  a  transformation  matrix 
which  outputs  the  LOS  coordinates,  “x”  and  “y”,  which  indicate  the 
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alignment/misalignment  of  the  secondary  mirror  assembly  with  the  primary  mirror 
assembly.  The  Structure  block  contains  the  structural  bending  dynamics.  As  shown  in 
Figure  3-3,  the  Structure  block  does  not  contain  the  augmented  rigid-body  dynamics.  The 
addition  of  rigid-body  states  will  be  shown  later  in  the  chapter.  The  MMAC  block 
represents  the  entire  moving-bank  multiple  model  algorithm  for  providing  active  control  of 
the  SPICE  structure. 


Figure  3-3.  SPICE  System  Block  Diagram  [2:  3-5] 


3.3.1  Disturbances 

The  structural  disturbances  are  modeled  as  entering  the  structure  at  six  points:  at 
the  top  and  bottom  of  each  tripod  leg.  Details  of  how  these  disturbances  impact  the 
structure  are  given  shortly  in  Equation  (3.3).  These  are  modeled  as  shaped  noise  over  the 
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31.416  to  62.832  radians  per  second  ( rps )  frequency  range.  This  is  modeled  as  6  white 
noise  inputs  of  equivalent  strength,  each  shaped  by  a  fourth-order  notch  filter,  as  shown  in 
Figure  3-4.  Thus,  the  disturbance  model  contributes  24  states  to  the  system  model.  The 
state  equation  for  this  model  is  given  as: 


XB(0  =  F„Xn(0  +  GnWB(f)  (3.1) 

where: 

•  XB  (t)  =  24-state  vector  representing  the  disturbance  states 

•  FB  =  24-by-24  constant  plant  matrix 

•  G„  =  24-by-6  constant  noise  input  matrix 

•  WB(f)  =  6-by-l  unit-strength  white  Gaussian  noise  vector 

The  corresponding  output  equation  is: 

n(f)  =  C„xB(f)  (3.2) 

where: 

•  n(0  =  6-by-l  output  colored  noise  vector 

•  CB  =  6-by-24  constant  matrix 


wr 


0.2905s2 

3947.8 

s2+ 44.422s +  985.95 

6 

s2+  88.84s +  3947.8 

Figure  3-4.  Disturbance  Model  Block  Diagram  [2:  3-6] 
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The  output  vector  is  defined  as: 


V 

"SMA  x  disturbance  Leg  1 

*2 

SMA  y  disturbance  Leg  2 

«3 

SMA  z  disturbance  Leg  3 

nA 

Bulkhead  z  disturbance  Leg  1 

Bulkhead  z  disturbance  Leg  2 

.n6_ 

Bulkhead  z  disturbance  Leg  3 

3.3.2  Proof  Mass  Actuators  (PMAs) 

Each  of  the  18  PMAs  are  modeled  identically  as  shown  in  Figure  3-5.  They  are 
modeled  as  second  order  high-pass  filters  with  a  cutoff  frequency  of  approximately  30  rps. 
Since  the  PMA  design  is  second  order,  the  18  PMAs  contribute  36  states  to  the  system 
model.  The  state  equation  for  the  PMA  model  is  given  as: 

*pma  (0  =  Fpmaxpma(0  +  BpmaU  pma(t)  +  GpmaVipma(t)  (3.4) 


where: 

•  xpma(0  =  3u-state  vector  representing  the  PMA  filter  states 

•  Fpma  =  36-by-36  constant  plant  matrix 

•  =  36-by-54  constant  deterministic  input  matrix 

•  Gpma  =  36-by-18  constant  noise  input  matrix 

•  Upma  --  54-by-l  “control”  vector  composed  of  components  Ucmd,  UvW  and 

•  W  pm(l)  =  18-by-l  unit-strength  white  Gaussian  noise  vector 
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Note  that  Ucmd  is  the  control  output  of  the  MMAC.  Additionally,  UvW  and  Upas  are  yve/ 
and  ypos,  respectively,  of  Figure  3-3,  i.e.,  outputs  of  the  “structure”  block.  These  are  the 

relative  velocities  and  positions  of  each  PMA  mass,  taken  with  respect  to  the  point  on  the 
structure  where  the  PMA  is  attached  [2:  Chapter  3  pg.  12]. 

The  corresponding  output  equation  is  given  by: 

V  (3.5) 
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where: 

•  ypBM(0  =  18-by-l  output  response  vector 

•  Cp^  =  18-by-36  constant  output  matrix 

•  Du  pma-  18-by-54  constant  deterministic  input  direct  feedthrough  matrix 

•  Dwpma  =  18-by-l 8  noise  direct  feedthrough  matrix 


The  direct  feedthrough  matrices  result  from  an  equal  order  over  equal  order  filter  transfer 
function.  The  control  and  noise  input  vectors,  upma  and  W  ,  respectively,  are  defined 

as: 


and 


'PMA  1  command 

PMA  2  command 

uis 

PMA  18  command 

U\9 

PMA  1  velocity 

U20 

= 

PMA  2  velocity 

^36 

PMA  18  velocity 

l/37 

PMA  1  position 

un 

PMA  1  position 

"54. 

PMA  18  position 

*  PMA  1  noise  " 

w  „„„  = 

pma 

w2 

= 

PMA  2  noise 

PMA  18  noise 

(3.6) 


(3.7) 
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The  output  vector,  yp(mj,  is  defined  by: 


*  PMA  1  response  " 

y2 

= 

PMA  2  response 

_y18_ 

PMA  18  response 

(3.8) 


3.3.3  Structure 

This  section  first  presents  the  equations  associated  with  the  flexible  structure.  This 
is  followed  by  the  augmentation  of  the  flexible  dynamics  with  the  rigid-body  dynamics. 
The  desired  coupling  is  formed  from  this  augmentation  in  order  to  excite  flexible  bending 
motion  of  the  structure  when  a  rigid-body  control  is  applied. 


3.3.3.1  Flexible  Dynamics 

In  general,  the  flexible  dynamics  of  the  structure  are  represented  by  the  following 
state  equation: 


xs(t)  =  Fsxs(t)+Bsypma(t)+Gsn(t)  (3.9) 

where 

•  Xs(t)  =  2n-by-l  state  vector  representing  the  structural  flexible  bending  modes 

•  n  =  number  of  structural  flexible  bending  modes 

•  Fs  =  2«-by-2n  matrix 

•  Bs  =  2n-by- 18  constant  matrix 
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•  Gs  =  2n-by-6  constant  matrix 

•  n  is  defined  by  Equations  (3.2)  and  (3.3) 

•  )fpma  IS  defined  by  Equations  (3.5)  and  (3.8) 

The  associated  output  equation  is  defined  as: 

y,(0  =  C,x,(0 + DnJr»(/) + Dyjy  ^(t) 


where 

•  ys  =  (f)  93-by-l  PMA  output  response  vector 

•  C5  =  93-by-2n  constant  output  matrix 

•  Dn5  =  93-by-6  noise  direct  feedthrough  matrix 

•  Dy  t  =  93-by-l 8  constant  deterministic  input  direct  feedthrough  matrix 


m  ji  jii  i  ...  I  9 


(3.10) 
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The  corresponding  output  vector  is  defined  as: 


LOS  optical  element  1,  yi<Ml 

y2 

LOS  optical  element  2,  ylos2 

y39 

LOS  optical  element  39,  y^39 

y40 

PMA  1  acceleration,  yafd 

y4i 

PMA  2  acceleration,  y^ 

y5? 

PMA  18  acceleration,  y^, g 

y58 

PMA  1  velocity,  yve/1 

y59 

PMA  2  velocity,  yvel2 

y74 

PMA  18  velocity,  yvein 

Yl5 

PMA  1  position,  y^, 

y7  6 

PMA  2  position,  ym2 
« 

y93. 

• 

PMA  18  position,  yposl 8 

The  first  39  elements  comprise  the  ytoj  vector,  as  indicated  in  Figure  3-3.  These  LOS 
outputs  are  obtained  from  an  optical  scoring  system  which  uses  39  separate  laser  and 
sensor  pairs  to  detect  changes  in  position  of  the  secondary  mirror  assembly  with  respect  to 
the  primary  mirror  assembly.  These  LOS  outputs  are  multiplied  by  a  LOS  transformation 
matrix  in  order  to  obtain  the  “x”  and  “y”  line-of-sight  errors,  as  indicated  in  Figure  3-3. 
The  yMC  output  vector  (elements  40  through  57)  indicates  the  relative  acceleration  of 

each  PMA  mass  with  respect  to  the  point  on  the  structure  where  the  PMA  is  attached. 
These  acceleration  outputs  are  used  by  the  accelerometer  model,  yet  to  be  discussed  in 
Section  3.3.4.  The  remaining  velocity  and  position  outputs  are  similar  to  the  acceleration 
outputs  and  were  discussed  in  Section  3.3.2. 
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As  indicated  by  Equation  (3.9),  the  flexible  structure  modeling  adds  2 n  states  to 


the  system  model.  The  next  section  will  address  the  augmentation  of  rigid-body  states. 


3.3.3.2  Rigid-body  Dynamics 

The  rigid-body  motion  considered  in  this  research  is  planar  motion.  Thus,  two 
additional  states  are  required.  In  polar  coordinate  terms,  these  are  angular  position  and 
angular  rate.  For  the  purposes  of  this  research  with  the  SPICE  2  flexible  structure,  this 
motion  is  confined  to  be  perpendicular  to  the  Z  axis,  rotating  about  an  axis  in  the  X-Y 
plane.  The  axis  of  rotation  is  located  at  the  base  (bulkhead)  of  the  structure.  Refer  to 
Figure  3-6  for  an  illustration  of  a  typical  rigid-body  rotation. 


Figure  3-6.  Rigid-Body  Rotational  Motion 
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The  two-dimensional  rigid-body  dynamics  can  be  described  by: 


mo = mo 


(3.12) 


where  A/(f)  is  the  moment  causing  the  rotation  to  occur,  /  is  the  constant  moment  of 
inertia  matrix  about  the  point  of  rotation,  and  0(f)  is  the  angular  acceleration.  Put  into 
state  space  form,  this  becomes: 
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where 

•  xXrb{t)  =  rigid-body  angle  position  =  0(f) 

•  x2 rb  ( 0  =  rigid-body  angular  rate  =  0(f) 

•  urb(0  =  rigid-body  control  applied  =  Af(f) 


In  more  general  terms.  Equation  (3.13)  can  be  written  as: 


(3.13) 


xrb(0  =  ^H,(0  +  hrburb(0 


(3.14) 


When  the  rigid-body  dynamics  are  augmented  to  the  flexible  dynamics,  Equation  (3.9) 
becomes: 


"*,(0  ' 

= 

X 

0 

'x,(f) ' 

+ 

X 

^  pma 

y  „.(«)' 

+ 

Gs 

_M0. 

0 

_*rb(0_ 

0 

.  Urb(< )  . 

0 

n(f) 


(3.15) 
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Note  that  an  additional  column  is  augmented  to  the  flexible  structure  B,.  The 
vector  is  simply  one  particular  column  of  the  B5  matrix.  This  comes  from  the  way  in 
which  the  rigid-body  control  is  applied  to  the  structure.  For  ease  of  simulation,  rigid-body 
control  is  applied  at  the  structure’s  base  at  one  of  the  PMA  locations.  (PMA  location 
number  13  was  used  in  this  research;  recall  Figure  3-1.)  In  order  for  the  flexible  dynamics 
of  the  system  to  be  realized  from  the  action  of  applying  a  force  to  a  particular  location  on 
the  structure,  the  appropriate  column  of  Bf  must  be  included  in  the  augmentation,  so  that 
it  is  correctly  multiplied  by  u^it). 


3.3.4  Accelerometers 

Each  of  the  18  accelerometers  is  modeled  as  a  first  order  high-pass  filter  with  an 
additive  corruptive  shaped  white  noise.  The  noise  filter  is  modeled  as  a  second  order 
high-pass  filter.  The  block  diagram  for  the  accelerometer  design  is  shown  in  Figure  3-7. 


Figure  3-7.  Accelerometer  Design 


With  18  inputs,  the  accelerometer  model  contributes  54  states  to  the  system  model.  Hie 
state  equation  for  the  accelerometer  design  is  given  as: 


X«cc(  0  = 


\(t)  ■ 

X  o ' 

'xa(t)  ' 

+ 

Ba 

Y^(0  + 

'  0  ' 

>na(0. 

.0  Fna_ 

_Xn  JO. 

0 

#  flfC'  ' 

_Gwa_ 

v»r  (f)  (3.16) 


where 

•  Xa(r)  =  18-state  vector  representing  the  accelerometer  response  states 

•  Xna(?)  =  36-state  vector  representing  time-correlated  accelerometer  noise 

•  Fa  =  18-by-18  constant  accelerometer  plant  matrix 

•  Fna  =  36-by-36  constant  accelerometer  noise  plant  matrix 

•  Bfl  =  18-by-18  constant  input  matrix 

•  Gwa  =  36-by-18  constant  noise  input  matrix 

•  Yaccit)  is  defined  by  elements  40  through  57  of  Equation  (3. 1 1) 

•  W 0^  =  18-by-l  unit-strength  white  Gaussian  noise  vector 


The  output  of  the  accelerometers  is  the  measurements  and  is  given  by  the  measurement 
equation: 

Z</,)  =  [H„  H.„] 

where 

•  Xa(t'  )=  18-state  vector  representing  the  accelerometer  response  states 

•  Xna  (ti )  =  36-state  vector  representing  time-correlated  accelerometer  noise 

•  Ha  =  18-by-l 8  constant  accelerometer  measurement  matrix 


*na(0. 


+Dyaccyacc(ti)+Dwawacc(ti) 


(3.17) 
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•  Hna  =  18-by-36  constant  accelerometer  noise  measurement  matrix 

•  Dya  =  18-by- 1 8  constant  direct  feedthrough  matrix 

•  Dwa  =  36-by-18  constant  direct  feedthrough  matrix 

•  YacM  is  defined  by  elements  40  through  57  of  Equation  (3.1 1) 

•  wocc  =  18-by- 1  unit-strength  white  Gaussian  noise  vector 

The  direct  feedthrough  matrices  in  the  second  and  third  terms  of  Equation  (3.17)  are  a 
result  of  the  accelerometer  model  containing  transfer  functions  with  equal  order  numerator 
and  denominator. 


3.4  Truth  Model  and  Filter  Model 

The  system  model  presented  thus  far  contains  114  +  2 n  states.  The  1 14  states  are 
comprised  of  PM  A,  accelerometer,  and  input  disturbance  and  accelerometer  measurement 
noise  states.  The  structure  contributes  the  2 n  states  for  n  flexible  modes.  The  truth  model 
used  by  Gustafson  [2]  will  be  employed  here,  and  that  is  with  a  40-mode  structure  model. 
This  gives  194  states  for  the  truth  model.  Using  40  flexible  bending  modes  as  the  truth 
model  creates  a  model  with  undamped  natural  frequencies  below  100  Hz,  which  is  the 
frequency  limit  of  the  research  model.  A  194-state  truth  model  is  also  computationally 
feasible  for  simulation  purposes. 

The  filter  model  consists  of  the  same  truth  model  states  with  the  shaping  noise 
states  all  replaced  by  white  noise.  This  reduces  the  number  of  states  to  54  +  2 n  =  134 
with  40  flexible  bending  modes.  Gustafson  also  investigated  other  reduced-order  filter 
model  designs  based  on  the  techniques  of  modal  reduction  and  internally  balanced  model 
reduction.  The  MMAC  did  not  perform  well  with  these  reduced-order  filter  designs, 
however.  This  poor  controller  performance  may  have  been  partially  due  to  the  low 
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number  of  bending  modes  utilized  (6  bending  modes).  Current  research  by  Schiller  [28] 
with  the  subsequent  SPICE  4  model  has  utilized  26  bending  modes  in  a  reduced-order 
filter  and  controller,  which  has  produced  excellent  MMAC  performance.  While  not  all  of 
the  improvement  in  controller  performance  may  be  attributed  to  the  higher  number  of 
bending  modes  utilized,  this  does  certainly  warrant  another  look  with  the  SPICE  2  model. 
Thus,  another  reduced-order  filter  model  consisting  of  26  bending  modes  (52  bending 
mode  states)  +  54  PMA  and  accelerometer  states  (106  total  states)  will  be  considered  in 
this  research.  The  reduced-order  filter  model  will  be  constructed  using  the  modal 
reduction  method  described  in  Chapter  2.  Particularly  once  more  precise  measurements 
are  assumed,  this  reduced-order  design  model  may  well  yield  good  performance. 


3.5  Summary 

This  chapter  has  presented  a  brief  discussion  of  the  SPICE  2  system  model  used 
for  this  research.  This  included  a  description  of  the  physical  structure,  along  with  a  look 
at  each  subsystem.  The  mathematical  models  for  input  disturbances,  PMAs,  the  flexible 
and  rigid-body  dynamics  and  accelerometer  measurements  were  presented.  The  truth 
model  was  shown  to  consist  of  194  states  while  the  filter  models  were  shown  to  consist  of 
134  and  106  states. 

The  next  chapter  will  discuss  all  aspects  of  the  various  software  programs  used  for 
simulating  the  MMAE  and  MMAC  algorithms. 
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4.  Simulation 


4.1  Introduction 

This  chapter  presents  a  discussion  of  the  software  and  the  computer  simulation 
plan  used  to  analyze  the  effectiveness  of  the  moving-bank  MMAE  and  MMAC  algorithms 
as  applied  to  the  SPICE  2  structure  model.  A  Monte  Carlo  analysis  is  the  main  instrument 
used  to  show  this  effectiveness  and  is  performed  using  the  MMAE/MMAC  software 
developed  from  previous  thesis  research  [1,  2,  5,  7,  1 1,  26,  29,  32].  The  formation  of  the 
discretized  parameter  space  and  the  changes  made  to  the  existing  software  will  also  be 
discussed  in  this  chapter. 


4.2  Monte  Carlo  Analysis 

A  Monte  Carlo  analysis  is  performed  in  order  to  evaluate  the  performance  of  the 
moving-bank  MMAE  and  MMAC  algorithms.  This  is  accomplished  by  using  the 
simulation  software  (see  Section  4.3)  to  produce  multiple-run  sample  statistics  of  the  data 
of  interest  (estimation  errors  for  filter  performance  and  line-of-sight  deviations  for 
controller  performance,  as  well  as  parameter  estimation  errors  for  both  cases).  The  main 
emphasis  of  the  simulations  will  be  with  the  194-state  truth  model  versus  the  134-state 
filter  model,  although  one  comparison  is  made  with  a  106-state  reduced-order  filter  model 
(see  Chapter  3). 

The  simulations  performed  can  be  broken  into  two  areas:  estimator  analysis  and 
controller  analysis.  These  two  are  depicted  in  Figure  4-1  with  variables  defined  as: 
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Figure  4-1.  (a)  Estimator  Simulation,  and  (b)  Controller  Simulation  [29] 


•  X,  (ti )  =  truth  model  states 

•  Xfitj )  =  filter  estimates  of  the  system  states 

•  a(  (t, )  =  parameter  vector  implemented  in  the  truth  model 

•  a^(t, )  =  filter  estimates  of  the  uncertain  parameter  vector 

•  eo  (h ) =  error  the  parameter  estimate  defined  as:  ea  (t, )  =  a,  {ti  )-a/(tl  ) 

•  ex(ti)  =  error  in  the  system  estimate 
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The  following  sections  present  the  error  vector  formulation  and  the  error  vector 


statistics. 


4.2.1  Error  Vector  Formulations 

The  X  and  Y  Line-of-Sight  (LOS)  deviations  are  the  concern  in  this  research. 
Therefore,  for  the  estimator  performance  evaluation,  an  error  vector  is  formed  by 
subtracting  the  filter-determined  X  and  Y  LOS  values  from  the  corresponding  truth- 
model-calculated  values.  This  is  shown  in  Figure  4- 1(a)  and  by  the  following  equation: 

K 

e,(»,)  =  C,x,(ri)-^C/ji//(,)pJ(ri)  (4.1) 

J=1 


where  C,  and  Cf.  are  the  output  matrices  used  to  determine  the  X  and  Y  LOS  deviations 
for  the  truth  model  and  /h  filter  model  respectively.  The  summation  in  this  equation  is 
computed  by  use  of  the  conditional  probabilities  from  each  of  the  filters  in  the  moving 
bank  (as  presented  in  Chapter  2).  The  actual  form  of  this  error  vector  is  given  by: 


X  -  axis  LOS  position  error 

Y  -  axis  LOS  position  error 

(4.2) 


For  the  controller  performance  evaluation  (closed  loop  performance),  the  item  of 
primary  concern  is  how  well  the  LOS  deviations  are  quelled  to  zero.  Thus,  the  vector  of 
importance  is  given  by  the  following  equation: 

e'x(ti)  =  C,X,(ti)  (4.3) 
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where  6 '(*,  )  is  the  vector  of  the  true  X  and  Y  LOS  deviations  (position  errors).  The 
control  vector  u(ti )  in  Figure  4- lb  also  warrants  investigation,  to  ensure  that  commanded 
controls  are  of  reasonable  magnitude.  In  both  estimation  and  control,  ea  of  Figure  4-1 
also  merits  inspection,  to  determine  how  effectively  the  algorithm  estimates  the  uncertain 
parameter  a  in  its  attempt  to  provide  adaptive  estimation  or  control  of  states. 


4.2.1  Error  Vector  Statistics 

The  mean  or  average  value  of  the  error  vectors  just  defined  is  one  of  two  critical 
statistics,  and  the  main  statistic  used  for  portraying  parameter  estimation  performance,  in 
this  research.  The  mean  is  given  by:  [19:  129]: 

1  L 

£{e,«,)}  -  (4.4) 

L  1= 1 


where  L  is  the  number  of  Monte  Carlo  runs  made,  and  ex/  (f,  )  is  the  value  of  the  error 
vector  during  the  Ith  simulation  at  run  time  tr  The  covariance  of  the  error  signal  is 
calculated  by  [19:  130]: 


(fi )  =  £{[M<, )  -  (<,-  )}][•,('/  )-£{*,  (</  )}]T  | 

1  L 

= X  ('■  >}  -  ITT  M».  >  ‘4-5> 


Equations  (4.4)  and  (4.5)  represent  general  formulas  for  calculating  error  vector 
statistics  and  are  used  to  calculate  the  mean  and  covariance  of  e'(f,  )  and  ea(f,  )  as  weU. 
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Additionally,  this  research  will  use  L  =  10  Monte  Carlo  runs  to  generate  the  sample 
statistics. 

The  error  vector  statistics  described  above  will  be  useful  for  creating  plots  of  the 
desired  variables.  However,  in  order  to  quantify  the  performance  of  the  controller  in  a 
compact  manner,  suitable  for  tabulating,  another  statistic  is  used.  This  is  the  temporal 
average  of  the  RMS  value  and  is  calculated  by  [2:  Chapter  4  pg.  4]: 

=  j-  £  +  for  p  =  1  and  2  (4.6) 

™  jm-N+l  PP 

where  is  the  temporal  average  of  the  pth  component  of  e'  and  N  is  the  number  of 
sample  periods  over  which  the  temporal  averaging  is  accomplished.  For  this  research,  the 
temporal  average  will  be  taken  over  the  last  three  seconds  of  the  run  duration. 


4.3  Simulation  Software 

In  the  most  recent  study  by  Gustafson  [2],  the  software  (composed  of  MATRIXx 
[16]  and  FORTRAN  programs)  which  was  “handed  down”  from  previous  research  [1,5, 
7,  1 1,  26,  29,  32]  was  modified  for  the  SPICE  2  model.  The  software  is  still  resident  on 
the  Sun  workstations  and  is  split  into  three  sections:  (1)  preprocessor,  (2)  processor,  and 
(3)  postprocessor.  Each  of  these  is  discussed  below  as  to  its  function  and  the  changes 
made  during  this  research. 


4.3.1  Preprocessor 

Inherited  Preprocessor.  The  preprocessor  is  actually  two  separate  processors:  a 
MATRIXx  portion  and  a  FORTRAN  portion.  MATRIXx  is  used  to  generate  the 
discretized  truth  and  filter  model  matrices  as  well  as  the  steady-state  filter  and  controller 
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gain  matrices.  These  matrices  are  then  stored  in  files  (one  file  for  each  truth  and  filter 
model  based  on  a  particular  discretized  parameter  value)  to  be  read  in  by  the  preprocessor 
FORTRAN  program,  which  merely  “collects”  all  the  information  for  each  parameter 
model  and  stores  the  variables  into  two  files  (one  for  the  truth  models  and  one  for  the 
filter  models)  to  be  read  by  the  processor. 

Current  Preprocessor.  During  this  research  it  was  found  that  the  MATRIXx 
routine  for  obtaining  the  steady-state  controller  matrix,  G*,  did  not  converge  to  the 

proper  value.  Investigation  into  this  discrepancy  showed  that  the  MATRIXx  routine  did 
converge  properly  for  the  simple  example  in  the  MATRIXx  manual.  Further  investigation 
showed  that  the  MATLAB  [15]  routine  for  solving  the  steady-state  gain  matrix  did 
converge  properly  for  the  SPICE  2  model.  Thus,  all  the  preprocessor  MATRIXx  coding 
was  converted  to  MATLAB.  As  will  be  shown  in  Chapter  5,  even  though  the  MATRIXx 
routine  failed  to  converge  to  the  proper  controller  gain,  identical  results  for  the  closed- 
loop  temporal  average  values  were  obtained  as  compared  to  Gustafson’s  [2]. 
Additionally,  the  residual  covariance  calculation  and  its  inverse,  Ak  and  A*1  of  Equation 

(2.29),  were  removed  from  the  FORTRAN  portion  of  the  code  and  put  into  MATLAB. 
Problems  persisted  in  the  past  with  this  matrix  inverse  calculation,  but  the  MATLAB 
routine  showed  no  problems.  One  other  major  change  to  the  preprocessor  was  to  convert 
the  FORTRAN  code  to  double  precision  arithmetic  to  allow  for  greater  numerical 
accuracy.  The  double  precision  routines  were  compiled  with  the  IMSL  library  [6]. 

4.3.2  Processor 

Inherited  Processor:  This  set  of  FORTRAN  subroutines  performs  the  Monte  Carlo 
simulation.  This  is  done  by  propagating  the  truth  and  filter  models,  with  dynamics  and 
measurement  noises  simulated  by  appropriate  scaling  of  the  outputs  of  a  unity  variance 
random  number  generator  available  in  FORTRAN.  The  different  types  of  bank  movement 
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and  MMAE/MMAC  options  are  performed  as  well.  The  various  variables  to  be  analyzed 
are  output  to  separate  files  for  each  sample  time  and  run.  These  files  are  then  read  in  by 
the  postprocessor. 


Current  Processor.  A  few  changes  were  made  to  the  processor  for  enhancement 
and  to  allow  for  the  analysis  required  in  this  research.  First,  the  conversion  to  double 
precision  was  made,  as  with  the  preprocessor.  Second,  code  was  added  for  applying 
open-loop  dither  inputs.  Options  were  made  available  via  an  input  data  file  for  selecting 
square  wave,  sine  wave  (fixed  and  swept  frequency),  or  wide  band  noise.  Third,  the 
program  was  modified  to  allow  for  rigid-body  states  to  be  included,  augmented  to  the 
bending  mode  states  originally  simulated.  This  was  accomplished  for  the  purpose  of 
simulating  rigid-body  slew  commands  in  order  to  analyze  the  effects  on  uncertain 
parameter  identification  enhancement. 

Two  other  changes  were  made  to  the  processor  due  to  inadequacies  discovered. 
The  first  concerns  the  bounds  placed  on  the  probability  calculation.  From  the  discussion 
in  Chapter  2,  there  must  be  a  lower  bound  placed  on  the  calculated  pk  probabilities.  The 

value  for  this  bound  was  set  to  0.1,  which  was  too  high.  The  lower  bound  was  reset  to 
0.01.  Additionally,  it  was  found  that  the  code  allowed  for  an  upper  bound  to  be  placed  on 
the  calculated  probability  as  well.  This  upper  bound  was  removed  all  together.  The 
%^pcond  change  concerns  the  probability  density  calculation.  It  was  found  that,  if  the 
residuals  were  too  large,  the  exponential  term  in  Equation  (2.27)  was  driven  to  zero,  due 
to  finite  word  length  with  the  digital  computer.  In  many  situations,  this  would  occur  for 
each  of  the  k  filters,  resulting  in  the  entire  denominator  of  Equation  (2.26)  being  evaluated 
as  zero,  preventing  the  probabilities  from  being  calculable.  So,  a  check  was  put  in  the 
code,  such  that,  if  the  density  calculation  produced  zero,  it  was  reset  to  a  very  small  value 
(to  approximate  a  near-zero  value  and  allow  the  probability  calculation  to  continue). 
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4.3.3  Postprocessor 

Inherited  Post  Processor.  This  FORTRAN  program  calculates  all  the  pertinent 
statistics  used  for  analyzing  the  MMAE  and  MMAC  performance.  These  statistics  are 
output  to  separate  Hies  to  be  used  by  a  compatible  plotting  routine. 

Current  Post  Processor.  The  only  change  to  the  post  processor  was  the  conversion 
to  double  precision.  MATLAB  was  used  to  plot  the  data  output  by  the  post  processor. 

4.4  Analysis  Plan 

This  research  effort  is  directed  towards  improving  the  parameter  identification 
performance  by  use  of  various  dither  input  signals  prior  to  application  of  closed-loop 
control.  The  secondary  issues  are  threefold:  (1)  possible  satisfactory  performance  with  a 
26-bending  mode  reduced-order  filter  model,  (2)  possible  controller  performance 
improvement  via  better  accelerometer  measurements  (less  measurement  noise),  and  (3) 
possible  estimator/controller/parameter  identification  performance  improvement  via  an 
alternate  discretized  parameter  space  (different  from  that  used  by  Gustafson  [2]). 

Before  any  of  these  issues  are  explored,  however,  the  software  and  system  models 
will  be  verified  by  attempting  to  reproduce  the  single  filter  benchmark  results  obtained  by 
Gustafson  [2].  Additionally,  throughout  the  dither  signal  analysis,  a  comparison  will  be 
made  between  the  ME/A  and  ME/I  methods  of  probability  density  calculation  (discussed 
in  Chapter  2).  The  next  three  sections  discuss  in  further  detail  the  issues  concerning  input 
dither  frequency,  improve®  accelerometer  measurements,  and  parameter  space 
discretization. 

4.4.1  Dither  Frequency  Selection 

Based  on  the  research  accomplished  by  Hanlon  [4],  the  flexible  structure  modal 
frequencies  are  the  basis  for  selecting  the  dither  frequencies.  The  modal  frequencies 


4-8 


correspond  to  the  system  eigenvalues  and  are  useful,  since  inputs  based  on  these 
frequencies  should  excite  the  system  more  than  by  using  some  other  arbitrary  frequency. 
The  modal  frequencies  for  the  first  30  modes  are  listed  in  Table  4-1.  In  order  to  determine 
which  frequency  or  frequencies  are  best  suited  for  exciting  the  system,  the  dither  input 
frequency  will  be  swept  through  the  range  of  modal  frequencies  using  a  sinusoid  signal.  It 
is  expected  that  a  dither  frequency  corresponding  to  one  of  the  lower  bending  modes  will 
excite  the  system  most  effectively,  due  to  the  physical  characteristics  of  the  space  structure 
-  it  being  two  rigid  bulkhead  assemblies  attached  by  tripod  legs.  Once  this  “best”  dither 
frequency  is  determined,  it  will  be  used  as  the  basis  for  square  wave  and  wide  band  dither 
inputs  as  well. 


4.4.2  Improved  Accelerometer  Measurements 

In  this  analysis,  an  attempt  will  be  made  to  provide  a  benchmark  with  a  “perfect” 
measurement  model.  This  will  be  accomplished  by  changing  the  gain  in  Figure  3-7  to 
zero.  This  will  have  the  afreet  of  removing  all  accelerometer  measurement  noise.  Next, 
the  original  gain  will  be  reduced  by  one  and  two  orders  of  magnitude  in  order  to  simulate 
a  possible  realistic  improvement  in  the  measurement  device.  In  each  of  these  cases  it  will 
be  necessary  to  modify  (retune)  the  original  matrix  values  for  optimal  results. 

From  Table  4-1,  note  that  the  first  26  modal  frequencies  are  less  than  50  Hz. 
Recall  from  Section  2.7.2  that,  with  a  sampling  time  of  0.01  sec,  the  maximum  useful 
sinusoid  dither  frequency  is  50  Hz.  Thus,  the  dither  input  frequencies  used  in  this  research 
will  be  limited  to  the  first  26  modal  frequencies.  This  may  not  be  an  actual  limitation, 
however,  since  it  is  expected  that  the  optimum  dither  frequency  will  reside  within  the 
lower  modal  frequencies. 
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Mode 

Eigenvalue 

Damping  Ratio 

Natural  Frequency  (Hz) 

1 

-0.09729  ±  48.64j 

0.002 

7.7418 

2 

-0.1115  ±55.77j 

0.002 

8.8762 

3 

-0.1115  ±55.77j 

0.002 

8.8763 

4 

-0.2318  ±115.9j 

0.002 

18.446 

5 

-0.2318  ±115.9j 

0.002 

18.449 

6 

-0.2350  ±117.5j 

0.002 

18.720 

mm 

-0.3357  ±  167.8j 

0.002 

26.712 

8 

-0.3604  ±  180.2j 

0.002 

28.681 

9 

-0.3605  ±  180.2j 

0.002 

28.684 

10 

-0.3793  ±  189.6j 

0.002 

30.180 

11 

-0.4103  ±  205.2j 

0.002 

32.654 

12 

-0.4106  ±205.3j 

0.002 

32.673 

13 

-0.4164  ±208.2j 

0.002 

33.135 

14 

-0.4166  ±208.3j 

0.002 

33.150 

15 

-0.4168  ±208.4j 

0.002 

33.168 

16 

-0.4624  ±  23 1.2j 

0.002 

36.798 

17 

-0.4626  ±  23 1.3j 

0.002 

36.816 

18 

-0.4689  ±  234.4j 

0.002 

37.312 

19 

-0.5407  ±270.3j 

0.002 

43.025 

20 

-0.5409  ±  270.4j 

0.002 

43.040 

21 

-0.5628  ±  28 1.4j 

0.002 

44.784 

22 

-0.5697  ±284.8j 

0.002 

45.334 

23 

-0.5698  ±  284.9j 

0.002 

45.343 

24 

-0.6116±305.8j 

0.002 

48.669 

25 

-0.6215  ±310.7j 

0.002 

49.455 

26 

-0.6216  ±  3 10.8j 

0.002 

49.467 

27 

-0.6724  ±  336.2j 

0.002 

53.505 

28 

-0.6894  ±  344.7j 

0.002 

54.864 

29 

-0.7004  ±  350.2j 

0.002 

55.732 

30 

-0.7008  ±  350.4j 

0.002 

55.770 

Table  4-1.  Modal  Eigenvalues  and  Natural  Frequencies  for  the  First  30  Modes 
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4.4.3  Parameter  Space  Discretization  (Sensitivity  Analysis) 

A  sensitivity  analysis  is  performed  in  order  to  produce  a  set  of  discrete  parameter 
points.  This  is  accomplished  by  setting  the  truth  and  filter  models  to  the  nominal 
parameter  value,  i.e.,  the  nominal  values  for  the  undamped  natural  frequencies  of  the 
bending  modes.  Simulations  are  then  run  in  this  single  filter  mode  with  closed  loop 
control  applied.  On  each  successive  run,  a  scalar  multiplier  on  all  the  undamped  natural 
frequencies  in  the  truth  model  is  varied  by  a  small  amount,  8 ,  until  the  closed  loop  control 
becomes  unstable.  The  change  in  each  undamped  natural  frequency  is  given  by: 

©^=d  +  S)®oW  (4-7) 

Once  the  value  for  8  is  found  which  produces  unstable  control,  a  new  parameter  point  is 
found  by  reducing  this  value  for  8  by  10  %.  Then,  the  truth  and  filter  models  are  reset 
with  the  new  values  for  the  undamped  natural  frequencies,  and  the  process  is  repeated  in 
order  to  obtain  a  new  value  for  8  and  a  new  parameter  point.  This  entire  process  is 
repeated  until  the  desired  range  in  undamped  natural  frequency  is  obtained.  Parameter 
points  are  also  generated  by  changing  the  modal  frequencies  in  the  opposite  direction  by: 

©,^=(l-8)©oW  (4.8) 

This  method  of  performing  a  sensitivity  analysis  to  obtain  the  discrete  parameter 
points  allows  for  the  widest  spread  between  each  point,  without  any  two  adjacent  points 
being  so  far  apart  so  as  to  produce  unstable  closed-loop  control.  Also,  this  would  yield 
residuals  in  filters  at  adjacent  parameter  points  that  would  have  sufficiently  different 
characteristics  that  the  pk  computation  of  Equations  (2.26)  to  (2.29)  would  not  be 

incapacitated  by  such  residuals  being  virtually  indistinguishable  from  one  another. 
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However,  for  reasons  discussed  in  Chapter  5,  this  method  was  used  only  to  obtain  the  first 
point,  displaced  upward  in  frequency  from  the  nominal  values.  The  resulting  8  was 
0.0025.  Twenty-nine  parameter  points  were  then  generated  in  approximately  the  same 
span  as  used  by  Gustafson  [2:  Chapter  5,  pg.  35],  each  with  a  scalar  multiplier  spacing  of 
0.0025. 


4.5  Summary 

This  chapter  has  presented  a  discussion  of  the  simulation  analysis  to  be  performed 
in  this  research.  The  Monte  Carlo  simulation  has  been  described,  along  with  the 
appropriate  error  vector  formulations.  Also  included  in  this  discussion  has  been  the 
analysis  plan  which  has  shown  the  main  points  of  the  research  to  be  accomplished,  along 
with  a  detailed  look  at  the  formation  of  the  discretized  parameter  space. 

The  next  chapter  will  discuss  the  results  of  this  research. 
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5.  Results 


5.1  Introduction 

Research  was  conducted  in  accordance  with  the  analysis  plan  presented  in  Chapter 
4.  This  research  was  based  on  the  theory  and  procedures  presented  in  the  first  three 
chapters  as  well.  This  chapter  presents  the  results  of  this  research.  As  discussed  in 
Chapters  1  and  4,  the  areas  to  be  investigated  are: 

•  Use  dithering  techniques  to  enhance  uncertain  parameter  identifiability,  prior  to 
engaging  closed-loop  control 

•  Investigate  the  effects  of  purposeful  rigid-body  slew  commands  in  regard  to 
enhancing  uncertain  parameter  identifiability 

•  Investigate  an  alternative  parameter  space  discretization 

•  Investigate  the  effects  of  increased  measurement  sensor  precision 

•  Investigate  the  effectiveness  of  a  106-state  modal  reduced-order  filter/controller 

Throughout  this  presentation  of  results,  reference  will  be  made  to  graphs  located  in 
Appendix  A.  These  figures  include  three  types:  (1)  is  the  mean  ±  one  standard  deviation 
(±a)  of  the  error  between  the  truth  model  values  and  filter  estimates.  This  is  computed 
by  Equations  (4.4)  and  (4.5),  utilizing  the  error  vector  in  Equation  (4.1),  and  is  referred  to 
as  estimation  error.  (2)  is  the  mean  ±  one  standard  deviation  (±ct)  of  the  actual  LOS 
error.  This  is  computed  by  Equations  (4.4)  and  (4.5),  utilizing  the  error  vector  in 
Equation  (4.3),  and  is  referred  to  as  LOS  error.  (3)  is  the  mean  of  the  filter-computed 
parameter  location.  The  true  parameter  variation  is  shown  coincidentally  so  as  to  indicate 
the  effectiveness  of  the  multiple  model  adaptive  algorithm  in  estimating  the  undamped 
natural  frequency  of  the  bending  modes. 
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5.2  MMAEIMMAC  “Best  Case ”  Performance 


This  section  presents  the  results  from  a  single  and  a  multiple  (three)  filter 
simulation  with  an  unchanging  parameter  value.  The  single  filter  represents  the  artificially 
informed  filter  (non- adaptive  case).  The  multiple  filter  unchanging  parameter  value  case  is 
presented  in  comparison,  to  demonstrate  the  ability  of  the  MMAE/MMAC  algorithm  to 
maintain  “lock”  on  a  good  estimate  of  parameters,  and  thus  to  maintain  state  estimation 
and/or  control  performance  that  is  comparable  to  the  artificially  informed  benchmark. 


5.2.1  Single  Filter  Baseline 

A  duplication  of  Gustafson’s  work  was  accomplished  with  the  194-state  truth 
model  vs.  the  134-state  filter  model.  The  purpose  of  this  simulation  with  a  single  filter 
was  two-fold.  First,  it  provides  a  verification  (in  part)  that  the  currently  modified  software 
is  performing  properly.  Second ,  it  provides  a  “best  case”  scenario  (the  filter  model 
parameters  exactly  represent  the  truth  model  parameters)  for  comparison  purposes.  The 
plots  of  Figure  A-l  show  the  estimation  error  and  the  LOS  error  for  the  X  and  Y  axes 
with  no  control  applied.  Similarly,  Figure  A-2  shows  the  same  error  plots  but  with  closed- 
loop  control  applied.  Note  that  control  is  applied  at  the  0.5  second  point.  This  was  done 
to  allow  good  state  (and  parameter)  estimates  to  be  achieved  (after  an  initial  transient) 
before  they  are  used  to  generate  feedback  control.  Clearly,  the  LOS  deviations  are 
quelled,  not  to  within  1  |i  radian,  but  to  a  significant  degree.  The  temporal  average  values 
of  the  X  and  Y  LOS  errors  are  24.8  ^.radians  and  26.1  jl  radians,  respectively,  with  no 
control,  and  4.38  p.  radians  and  4.88  jl  radians,  respectively,  with  control  applied.  The 
LOS  deviations  with  control  applied  are  precisely  the  same  as  that  obtained  by  Gustafson 
[2:  Chapter  5,  pg.  5-30]. 
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5.2.2  Multiple  Model  Baseline 


The  7-point  discretized  parameter  space  developed  by  Gustafson  [2]  is  used  here 
with  the  true  parameter  set  to  the  nominal  (center)  position.  A  bank  of  three  filters  is 
used,  initially  centered  on  the  true  parameter  position.  For  the  ME/A  case.  Figure  A-3 
shows  the  estimation  and  LOS  error  for  the  X  and  Y  axes  with  control  applied.  The 
temporal  average  values  of  the  X  and  Y  LOS  errors  are  4.S8  p  radians  and  5.10  p  radians, 
respectively.  For  the  ME/I  case,  Figure  A-4  shows  the  estimation  and  LOS  error  for  the 
X  and  Y  axes  with  closed-loop  control  applied.  The  temporal  average  values  of  the  X  and 
Y  LOS  errors  are  5.36  p  radians  and  6.15  pradians,  respectively.  It  is  seen  here  that  the 
ME/A  method  performed  with  15  %  and  17  %  better  control  on  the  X  and  Y  axes, 
respectively.  A  summary  of  the  MMAC  performance  is  presented  in  Table  5-1. 


Configuration 

Open-loop 

24.8 

26.1 

Single  Filter 

4.38 

4.88 

Multiple  Model  (ME/A) 

4.58 

5.10 

Multiple  Model  (ME/I) 

5.36 

6.15 

Table  5-1.  MMAC  LOS  Temporal  Average  RMS  Deviations 


5.3  Parameter  Estimation  Performance 

This  section  presents  the  results  from  utilizing  the  7-point  parameter  space 
developed  by  Gustafson  [2].  Throughout  the  discussion  in  this  section,  the  analysis  is 
conducted  without  closed-loop  control  applied.  The  only  control  applied  is  the  dither 
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signal  or  the  rigid-body  slew  command  (open- loop  control).  In  this  research,  the  dither 
signals  were  input  via  the  control  variable,  u(/, )  at  the  six  PMA  points  on  the  bulkhead  of 

the  structure  (locations  13  - 18).  This  seemed  to  be  a  logical  choice  for  dither  input,  given 
the  shape  of  the  SPICE  structure.  With  the  larger  mass  located  at  the  bulkhead,  any 
disturbances  entering  the  structure  from  the  base  would  cause  the  tripod  legs  to  be  more 
easily  shaken. 

5.3.1  BaseUne  (No  Dither) 

Before  evaluating  the  model  with  dither  signals,  part  of  Gustafson’s  work  was 
again  duplicated  in  order  to  establish  a  baseline  to  compare  against.  Here,  two  sets  of 
simulations  were  run.  In  the  first  set,  the  true  parameter  is  allowed  to  “jump”  to  a  new 
value  midway  in  the  simulation.  Many  different  configurations  were  tested.  These 
included  having  the  bank  centered  on  the  true  parameter  and  centered  at  different 
locations  in  the  parameter  space  at  the  start  of  the  run.  Figure  A-5  is  representative  of 
these  simulations  with  the  ME/A  method,  which  closely  compares  to  Gustafson  [2: 
Chapter  5,  pg.  5-48].  Figure  A-6  is  for  the  ME/I  case.  Both  of  these  methods  provide 
poor  parameter  tracking  performance,  but  the  ME/A  method  displays  better  performance 
than  the  ME/I  method.  The  second  set  of  simulations  consisted  of  allowing  the  true 
parameter  to  start  at  the  1st  position  and  increase  by  one  point  every  three  seconds,  up  to 
the  7th  position.  Again,  the  bank  of  filters  is  centered  at  various  starting  positions  on 
different  runs.  With  the  true  parameter  value  changing  in  this  fashion,  the  moving  bank 
was  unable  to  track  the  varying  parameter  well,  under  asy  bank  center  starting  position. 
Figure  A-7  is  representative  of  these  simulations  with  ME/A.  Figure  A-8  is  the  same 
simulation  but  with  the  ME/I  method.  Notice  that  neither  method  tracks  the  true 
parameter  very  well,  but  that  the  ME/A  method  does  somewhat  better  than  the  ME/I. 
This  is  consistent  with  the  better  MMAC  closed-loop  performance  with  the  MEA 
method. 
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5.3.2  Dither  Signal  Enhancement  Performance 

Swept  Sine  Wave :  Chapters  2  and  4  discussed  the  idea  of  a  “best”  or 
identification-enhancing  dither  frequency.  In  order  to  determine  this  frequency,  a  sine 
wave  dither  was  applied  with  the  dither  frequency  sweeping  through  a  range  of  values 
(swept  frequency  sine  wave).  The  frequencies  were  in  the  range  of  7  Hz  to  40  Hz.  Recall 
from  Chapter  4  that  this  is  the  range  of  the  lower  bending  mode  frequencies.  By  using  this 
swept  sine  wave  technique,  the  best  dither  frequency  was  determined  to  be  32  Hz,  which 
corresponds  approximately  to  the  11th  modal  frequency  (see  Table  4-1).  Figure  A-9 
indicates  the  moving  bank  “zeroing  in”  on  the  true  parameter  value  of  6.  The  simulation 
shown  starts  with  a  sine  wave  frequency  of  29  Hz  and  is  incremented  by  1  Hz  every  2 
seconds.  So,  when  a  frequency  of  32  Hz  is  reached  at  6  seconds,  the  algorithm  has 
“latched”  on  to  the  true  parameter  point  value  of  6.  The  same  result  of  32  Hz  was 
obtained  using  both  ME/A  and  ME/I  methods.  Since  this  identification-enhancing 
frequency  did  not  correspond  to  precisely  one  of  the  modal  frequencies,  many  different 
simulations  were  run  in  order  to  eliminate  the  possibility  of  there  being  a  better  dither 
frequency  than  32  Hz.  From  Figure  A-9,  it  would  appear  that  33  Hz  would  work  as  well, 
and  in  some  situations  this  was  true.  However,  a  frequency  of  32  Hz  performed  better 
most  of  the  time  with  different  signal  amplitudes  and  with  different  parameter  positions 
declared  as  the  true  value.  The  amplitude  best  suited  for  this  model  was  found  to  be  100 
Newtons  (N).  Small  deviations  from  this  value  did  not  affect  the  results,  but  too  small  an 
amplitude  prevented  proper  parameter  identification,  and  too  large  a  value  disturbed  the 
system  to  such  a  degree  that  prevented  proper  parameter  identification  as  well. 

Sine  Wave :  With  the  dither  frequency  of  32  Hz  determined,  the  full  ranging  true 
parameter  variation  simulation  discussed  in  Section  5.3.1  was  accomplished  with  a  sine 
wave  dither  input  signal.  Figures  A-10  and  A-ll  show  the  result  for  the  ME/A  and  ME/I 
methods,  respectively.  It  can  be  seen  that  the  ME/I  method  performs  better  than  the 
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ME/A  method.  This  is  in  contrast  to  the  better  ME/A  performance  with  no  dither  applied. 
Note  that  in  Figures  A-10  and  A- 11,  the  parameter  estimate  is  at  times  “locked”  cm  to  a 
value  which  is  one  off  from  the  true  parameter  value.  For  example,  in  Figure  A- 11,  for  the 
first  3  seconds,  the  parameter  estimate  is  indicating  a  value  of  0.5,  while  the  true 
parameter  value  is  1.  The  reason  for  this  “locking”  on  to  the  incorrect  parameter  value 
may  be  due  to  the  fact  that  the  “best”  dither-enhancing  frequency  is  slightly  different  for 
each  declared  truth  model  (based  on  a  different  parameter  value).  Remember  that  the 
dither  frequency  of  32  Hz  was  chosen  because  it  worked  well  for  the  entire  range  of 
parameter  values. 

Square  Wave:  The  results  of  the  square  wave  dither  input  simulations  are  shown 
in  Figures  A- 12  and  A- 13  for  the  ME/A  and  ME/I  methods,  respectively.  As  with  the  sine 
wave  dither  signal,  the  ME/I  method  outperforms  the  ME/A  method  in  parameter 
tracking.  Many  simulations  were  run  with  different  square  wave  frequencies  (similar  to 
the  swept  frequency  sine  wave  simulations).  Still,  the  32  Hz  prevailed  as  the  parameter¬ 
enhancing  frequency.  It  was  seen  that,  even  with  the  extra  high  frequency  content,  the 
square  wave  did  not  adequately  enhance  parameter  identification,  when  other  than  a  32  Hz 
frequency  was  used.  Furthermore,  by  comparing  Figure  A-10  with  A- 12,  and  Figure  A-l  1 
with  A-13,  it  was  seen  that  there  was  no  substantial  difference  between  the  tracking 
performance  with  sine  wave  or  square  wave  dither  inputs. 

Wide  Band  Noise :  Frequencies  in  the  range  of  7  Hz  to  40  Hz  were  used  for  the 
wide  band  noise  input  signal.  However,  performance  was  degraded  when  frequencies 
above  33  Hz  were  used.  Different  values  for  the  gain  were  evaluated  (see  Figure  2-4) 
with  a  value  of  100  giving  the  best  performance.  Figure  A-14  indicates  the  tracking 
performance  of  the  ME/I  method  with  a  signal  frequency  range  of  7  Hz  to  33  Hz.  This 
frequency  range  was  then  changed  to  30  Hz  to  33  Hz.  Figures  A-15  and  A-16  show  the 
parameter  tracking  performance  of  the  ME/A  and  ME/I  methods,  respectively,  for  this 
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frequency  range.  Again,  as  with  the  sinusoidal  dither  input  signals,  the  ME^I  method 
provided  bettor  tracking  performance  than  the  ME/A  method. 


5.3.3  Rigid-body  Motion  Enhancement  Performance 

Before  the  rigid-body  matrices  could  be  augmented  to  the  modal  matrices  (see 
Section  3. 3. 3. 2),  a  moment  of  inertia  value  had  to  be  determined  for  use  in  Equation 
(3.13).  From  the  physical  description  of  the  SPICE  2  model  in  the  SPICE  Preliminary 
Design  Review  (PDR)  [14],  the  moment  of  inertia  was  determined  to  be  10000  kgm2. 
The  simulation  proceeded  with  a  rigid-body  command  input  of  +1000  N  for  1  second 
followed  by  an  input  of  -1000  N  for  1  second.  These  command  inputs  simulate  the 
rotating  structure,  accelerating  and  then  decelerating,  producing  a  total  slew  angle  of 
approximately  3.5  degrees.  This  is  consistent  with  the  test  performed  on  the  SPICE 
model  as  indicated  in  the  SPICE  PDR  [14:  TV-19]  and  thus  helps  to  verify  the  moment  of 
inertia  calculation.  Figures  A- 17  and  A- 18  show  the  parameter  estimation  performance 
with  the  slew  commands  just  described  for  the  ME/A  and  ME/I  methods,  respectively.  As 
these  figures  indicate,  the  rigid-body  motion  causes  a  significant  disturbance,  and  as  the 
disturbance  quells,  the  MMAE  is  able  to  track  the  true  varying  parameter  value, 
somewhat.  The  large  disturbance  shown  in  these  figures  occurs  between  the  1.5  second 
point  and  the  3.5  second  point  (the  time  period  the  slew  commands  were  applied).  Once 
the  rigid-body  slewing  motion  has  stopped,  the  algorithm  takes  approximately  0.5  seconds 
for  acquisition  of  a  good  parameter  estimate.  Note  that  the  excellent  parameter 
identification  performance  found  with  purposefully  constructed  dither  signals  is  not  found 
here,  but  that  the  parameter  identification  is  enhanced,  as  compared  to  the  performance 
with  no  dither  input  at  all.  Also  note  that  the  ME/A  and  ME/I  methods  perform  about  the 
same,  unlike  with  the  purposeful  excitation  forms  of  dither  enhancement 
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5.4  New  Parameter  Space  Diecretteatian 


The  simulations  discussed  in  Section  5.3  were  based  cm  the  7-point  parameter 
space  developed  by  Gustafson  [2].  However,  due  to  the  poor  performance  of  the 
parameter  estimate  tracking  without  dither  enhancement,  the  parameter  point 
discretization  was  re-investigated,  as  discussed  in  Section  4.4.3.  As  indicated  in  Chapter 
4,  the  first  parameter  point  obtained  was  only  a  factor  of  0.0025  from  the  nominal  point. 
Since  it  was  desired  to  maintain  the  same  range  of  undamped  natural  frequency  factors  as 
that  used  by  Gustafson,  this  new  parameter  discretization  would  require  many  more  than 
seven  points.  The  original  range  of  factors  was  0.9811  to  1.0466.  In  order  to  conserve 
time,  the  new  parameter  space  was  determined  by  equally  spacing  each  point  by  0.0025. 
The  new  range  used  was  0.98  to  1.05.  This  is  not  the  optimum  parameter  point  spacing, 
but  it  is  one  with  many  more  points  than  previously  used,  and  it  is  easily  generated.  For 
the  comparison  purposes  of  this  thesis,  this  will  suffice.  Thus,  the  new  discretized 
parameter  space  contains  29  points  (co,  through  OO^)  with  (o9  as  the  nominal  point,  as 

indicated  by  Table  5-2. 


Parameter 

©i 

®2  . 

■  (Og  m  • 

.  <oN 

Factor 

0.98 

0.9825  .  . 

.  1.0  .  . 

.  1.05 

Table  5-2.  Parameter  Space  Discretization 


The  same  set  of  simulations  as  presented  in  Section  5.3  were  run  with  the  new 
discretized  parameter  space  to  see  the  impact  of  the  fins'  discretization.  These  results  are 
presented  in  the  next  subsections. 
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5.4.1  Baseline  (No  Diiker) 


With  these  simulations,  the  true  parameter  varied  from  the  1st  to  the  29th  position, 
changing  position  every  second.  The  bank  of  filters  is  initially  centered  at  the  9th 
(nominal)  position.  Figures  A-19  and  A-20  represent  the  baseline  simulation  for  the 
ME/A  and  ME/I  methods,  respectively.  Note  the  greatly  improved  tracking  performance 
with  the  MDEVA  method  versus  the  ME/I  method.  This  result  is,  however,  consistent  with 
the  result  with  the  old  parameter  space  tracking  performance  in  that  the  ME/A  method 
performed  better  with  no  dither  enhancement 

5.4.2  Dither  Signal  Enhancement  Performance 

Sine  Wave :  Figures  A-21  and  A-22  indicate  the  parameter  tracking  performance 
of  the  ME/A  and  ME/I  methods,  respectively,  with  a  32  Hz,  100  N  amplitude  sine  wave 
dither  signal.  Notice  here  that  the  performance  rating  for  the  two  methods  is  mixed.  The 
ME/A  provides  better  tracking  in  the  middle  of  the  simulation  between  the  6th  and  24 
parameter  points,  while  the  performance  of  the  ME/I  method  is  better  at  both  ends. 

Square  Wave :  Figures  A-23  and  A-24  show  the  tracking  performance  of  the 
ME/A  and  ME/I  methods,  respectively,  with  a  square  wave  dither  input.  The  mixed 
performance  between  the  two  methods  is  very  similar  to  the  performance  with  the  sine 
wave  dither.  The  results  of  the  sine  wave  and  square  wave  dither-enhanced  performance 
are  also  consistent  with  the  results  from  the  7-point  parameter  space,  in  that  there  is  no 
substantial  difference  between  die  performance  with  sine  or  square  wave  dither  inputs. 

Wide  Band  Noise :  Figures  A-23  and  A-26  indicate  the  parameter  tracking 
performance  of  the  ME/A  and  ME/I  methods,  respectively,  with  wide  band  noise 
frequency  inputs  of  30  Hz  to  33  Hz.  A  gain  value  of  100  is  used,  as  with  the  old 
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parameter  space  simulations.  Again,  aa  with  the  sinusoidal  dither  inputs,  die  performance 
of  the  MEVA  and  MEVI  methods  is  mixed. 

Rigid-body  Motion :  The  same  rigid-body  commands  were  simulated  as  before. 
Figures  A-27  and  A-28  show  the  parameter  tracking  performance  of  the  ME/A  and  ME/I 
methods,  respectively.  Note  that  the  tracking  performance  of  both  these  methods  is  worse 
than  with  no  dither  input;  the  reason  for  this  phenomenon  is  not  clear. 


5.5  Summary  of  ME! A  vs  ME/I  Performance 

Table  5-3  indicates  the  compilation  of  the  comparison  study  accomplished  with  die 
ME/A  versus  the  ME/I  methods  of  probability  density  computation.  As  was  shown  in  die 


Performance  Measure 

Configuration 

ME/A 

Multiple  Model  MMAC 

Unchanging  Parameter 

mm 

Parameter  Tracking 

(7-point  Space) 

No  Dither 

V 

Sine  Wave 

mm 

Square  Wave 

mm 

Wide  Band  Noise 

■H 

Rigid-body  Motion 

mixed 

mixed 

Parameter  Tracking 
(29-point  Space) 

No  Dither 

V 

Sine  Wave 

■■ 

Square  Wave 

mm 

Wide  Band  Noise 

mixed 

mixed 

Rigid-body  Motion 

mm 

Table  5-3.  ME/A  versus  ME/I  Performance  (V  corresponds  to  better  performance) 
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preceding  discussion,  neither  method  is  better  in  every  situation.  In  Table  5-3,  a  “>/' 
(check)  indicates  which  method  was  better  in  die  given  configuration. 


5.6  Increased  Measurement  Precision  Performance 

Referring  to  Figure  3-7,  the  gain  of  1.47e-4  is  an  indication  of  the  variance  of  the 
truth  model  measurement  noise.  For  the  purpose  of  establishing  a  comparison  baseline, 
this  gal"  was  set  to  zero  (no  measurement  noise),  and  the  single  filter  closed-loop  control 
simulation  of  Section  5.2.1  was  run.  The  RMS  values  for  the  X  and  Y  LOS  deviations 
were  exactly  the  same  as  with  the  measurement  noise  included  (4.38  and  4.88  |i  radians, 
respectively).  Appropriate  re- tuning  of  the  Kalman  filter  R  matrix  values  was  attempted, 
with  no  change  in  the  result.  In  fact,  the  R  matrix  values  could  not  be  reduced  by  any 
significant  amount  due  to  the  [H(t,  )P(t,~  )HT  {ti )  4-  R(f, )]  term  in  Equation  (2.19) 
becoming  “too  small”,  preventing  the  steady  state  Kalman  filter  gain  from  being 
computable  due  to  difficulty  in  inverting  the  [H(f,  )P(^“  )HT  (/, )  +  R(f, )]  matrix.  Thus, 
increased  measurement  precision  had  no  effect  on  the  MMAC  performance.  Clearly,  more 
information  in  the  form  of  additional  measurements,  and/or  additional  actuators,  must  be 
included  before  desired  controller  performance  will  be  realized. 


5. 7  Reduced-Order  Design  Performance 

The  reduced-order  design  model  was  constructed  according  to  the  modal 
reduction  technique  described  in  Chapter  2.  The  modal  reduced  design  consisted  of  the 
first  26  bending  modes  (52  bending  mode  states)  plus  54  PMA  and  accelerometer  states, 
for  a  total  of  106  states.  This  106-state  reduced-order  filter  was  run  against  the  194-state 
truth  model  in  the  single  filter  configuration  with  closed-loop  control  applied.  Figure  A- 
29  shows  the  estimation  and  LOS  error  plots  with  closed-loop  control  applied.  Despite 
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the  relatively  high  number  of  bending  mode  states  included  in  this  filter  design,  the 
algorithm  was  unable  to  provide  sdequate  control  to  quell  the  vibrating  structure.  The 
temporal  averages  for  the  X  and  Y  LOS  errors  were  21.0  and  23.3  p  radians,  respectively. 
The  result  is  the  same  as  that  obtained  by  Gustafson  [2:  Chapter  S,  pg.  5-21]  with  his  6 
bending  mode  (66-state)  reduced- order  filter  model. 


5.8  Summary 

This  chapter  has  presented  the  results  from  the  analysis  conducted  in  this  research. 
First  presented  was  the  artificially  informed  single  filter  performance  baseline,  along  with  a 
verification  of  Gustafson’s  [2]  single- filter  result.  For  comparison  purposes,  the  multiple 
model  adaptive  algorithm  performance  with  an  unchanging  parameter  value  was  also 
shown.  The  results  of  implementing  the  various  dithering  techniques  were  then  presented. 
The  results  indicated  the  excellent  parameter  enhancement  with  the  sine  and  square  wave 
dither  signals.  The  wide  band  noise  dither  signal  was  also  found  to  provide  effective 
enhanced  parameter  identification.  The  results  of  naturally  exciting  the  structure  through 
the  use  of  rigid-body  slew  commands  did  improve  the  performance  of  parameter 
identification,  although  not  to  the  degree  exhibited  with  the  purposeful  dither  excitation. 
This  was  followed  by  a  presentation  of  results  of  parameter  identification  utilising  the  new 
discretized  parameter  space.  The  striking  difference  in  performance  with  the  new,  more 
finely  discretized  parameter  space  was  the  very  good  parameter  tracking,  with  no  dither 
enhancement.  Improved  parameter  identification  was  also  obtained  with  dither 
enhancement.  Throughout  the  discussion  of  MMAE/MMAC  and  parameter  estimation 
tracking  performance,  a  comparison  was  made  of  the  ME/A  and  ME/I  methods  of 
hypothesis  conditional  probability  computation.  As  was  shown  in  Table  5-3,  neither 
method  is  better  than  the  other  in  every  situation.  Lastly,  the  results  of  increased 
measurement  precision  and  reduced-order  filter  design  were  shown.  Although  these  last 
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two  areas  of  investigating  did  not  show  the  “positive”  results  as  that  obtained  with  the 
other  areas  of  this  thesis  research,  die  results  obtained  did  shed  light  into  the  limitations  of 
the  SPICE  2  model.  Chapter  6  will  discuss  the  conclusions  drawn  from  analyzing  the 
results  displayed  in  this  chapter,  as  well  put  forth  some  recommendations  for  further 
research  in  this  area. 
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6.  Conclusions  and  Recommendations 


6.1  Introduction 

The  main  purpose  of  this  thesis  has  been  to  research  and  design  effective  dither 
inputs  for  enhanced  parameter  identifiability  for  the  SPICE  2  flexible  space  structure. 
Previous  work  by  Gustafson  [2]  indicated  the  need  for  enhanced  parameter  identification. 
Thus,  this  research  with  the  SPICE  2  structure,  utilizing  moving-bank  Multiple  Model 
Adaptive  Estimation  and  Control  (MMAE/MMAC)  algorithms,  has  been  a  direct 
continuation  of  Gustafson’s  work.  Purposeful  excitation  dither  inputs  consisting  of 
sinusoidal  waves,  square  waves,  and  wide  band  noise  were  investigated,  as  well  as 
purposeful  rigid-body  slew  commands.  The  rigid-body  motion  was  applied  to  determine 
the  effects  on  enhanced  parameter  identification  through  natural  excitation  of  the 
structure’s  bending  modes.  Other  issues  explored  in  this  research  were:  (1)  a  new 
discretized  parameter  space,  (2)  increased  measurement  precision,  and  (3)  a  reduced-order 
filter/controller. 

The  ultimate  conclusions  reached  in  the  course  of  this  research  are:  (1)  dither- 
enhanced  parameter  estimation  performs  very  well  for  the  SPICE  2  flexible  structure,  and 
(2)  better  measurement  precision  and  reduced-order  filter/controller  designs  (up  to  26 
structure  modes)  do  not  yield  adequate  state  estimation  and  control  performance  for  the 
MMAE/MMAC  algorithms. 
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6.2  Conclusions 

The  analysis  conducted  in  this  research  began  with  a  verification  of  Gustafson’s 
work  with  the  single  filter/controller  design.  This  was  followed  by  establishing  baseline 
performances  with  the  7-point  parameter  space,  consisting  of  performance  against  (1)  an 
unchanging  parameter  value  and  (2)  a  varying  parameter  value,  both  with  no  dither 
enhancement.  Here,  it  was  seen  that  the  ME/A  method  performed  better  than  the  ME/I 
method  of  hypothesis  conditional  probability  computation.  In  contrast,  all  of  the 
performance  evaluations  conducted  with  purposefully  selected  dither  inputs  indicated  that 
the  ME/I  method  provided  better  parameter  identification  than  the  ME/A  method.  Of  the 
various  forms  of  dither  inputs  evaluated  for  parameter  identification  performance,  each 
performed  very  well.  While  the  sinusoid  dither  signals  were  better  able  to  “latch”  on  to 
the  true  parameter  value,  the  best  parameter-identification  dither  frequency  had  to  be 
correctly  determined.  With  the  wide  band  noise  dither  input,  the  “latching”  on  effect  was 
not  seen,  however,  a  broad  range  of  dither  frequencies  could  be  used  with  very  good 
performance. 

Thus,  each  type  of  dither  input  has  its  advantage,  and  a  combination  of  the  two 
types  could  be  used  in  a  real  application.  For  example,  wide  band  noise  dither  could  help 
identify  the  correct  range  of  frequencies.  The  noise  band  could  be  varied  to  the  point  that 
the  parameter-identification  frequency  would  be  known.  Then,  a  switch  to  a  sinusoid 
dither  at  this  frequency  could  be  made. 

Parameter  identifiability  enhancement  was  also  observed  from  the  effects  of  rigid- 
body  rotation-induced  excitation  of  the  bending  modes.  However,  the  effects  from  this 
type  of  commanded  slew  input  were  not  as  pronounced  as  those  from  the  sinusoid  and 
wide  band  noise  dither  inputs.  Indeed,  parameter  tracking  was  improved  with  rigid-body 
slewing  and  could  possibly  be  exploited  in  a  real  situation  whenever  the  structure  is  moved 
(pointed)  in  a  different  direction. 
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The  conclusion  drawn  from  simulations  conducted  with  the  new  29-point 
discretized  parameter  space  is  that  this  new  parameter  space  is  a  viable  alternative  to  the 
7-point  space,  if  no  dither  enhancement  is  to  be  used.  With  dither -enhanced  parameter 
identification,  the  two  parameter  space  discretizations  seemed  to  perform  equally  well. 
Thus,  dither  enhancement  could  be  used  with  the  7-point  parameter  space  as  an  alternative 
to  the  larger  memory  storage  requirement  with  the  29-point  parameter  space. 

The  analyses  performed  with  the  increased  measurement  precision  and  the 
reduced-order  filter/controller  revealed  some  limitations  with  the  SPICE  2  model.  The 
conclusions  drawn  from  these  results  are:  (1)  more  and/or  different  types  of  measurements 
may  be  necessary  for  improved  MMAC  performance,  (2)  better  and/or  additional 
actuators  should  be  considered. 


6.3  Recommendations 

Dither-enhanced  parameter  identification  was  shown  by  this  research  to  be  very 
effective  with  the  moving-bank  MMAE/MMAC  algorithms,  as  applied  to  the  SPICE  2 
flexible  structure.  Recommendations  for  further  research  are: 

1 .  Implement  the  latest  SPICE  structure  model. 

2.  Apply  the  dither-enhancing  methods  to  other  multiple  model  adaptive  research. 

3.  Investigate  the  poor  parameter  tracking  performance  of  the  29-point  parameter 
space  with  rigid-body  rotation-induced  excitation  of  the  bending  modes. 

4.  Consider  more  than  just  a  single  scalar  multiplier  on  bending  mode  frequencies 
as  the  uncertain  parameter.  Use  a  vector  of  parameters,  thereby  allowing  for 
different  relative  spacing  of  the  modal  frequencies  as  well  as  their  simultaneous 
movement  on  the  frequency  axis. 
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Appendix  A:  MMAE,  MM  AC,  and  Parameter  Estimation  Performance  Plots 

This  appendix  contains  all  the  plots  corresponding  to  the  results  discussion  in 
Chapter  5.  Figures  A-l  through  A-4  and  A-29  contain  mean  ±  one  standard  deviation 
(±o)  plots  of  estimation  error  and  LOS  error.  Each  of  these  four  figures  contain  four 
plots,  labeled  (a),  (b),  (c),  and  (d).  (a)  and  (b)  correspond  to  the  X-  and  Y-axis  estimation 
error  plots;  (c)  and  (d)  correspond  to  the  X-  and  Y-axis  LOS  error  plots.  The  remaining 
figures  are  plots  of  the  mean  parameter  estimate.  Shown  also  on  these  plots  is  the  true 
parameter  value,  for  comparison. 
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Figure  A- 1.  Single  Filter  Estimation  and  Control  Baseline  (Open  Loop) 
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Figure  A-2.  Single  Filter  Estimation  and  Control  Baseline  (Closed  Loop) 
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Figure  A-3.  Multiple  Filter  Estimation  and  Control  Baseline  (Closed  Loop,  ME/A) 
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Figure  A-4.  Multiple  Filter  Estimation  and  Control  Baseline  (Closed  Loop,  MEVI) 
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Figure  A-5.  7-Point  Parameter  Space:  True  Parameter  Jump,  No  Dither  (ME/A) 
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Figure  A-6.  7-Point  Parameter  Space:  True  Parameter  Jump,  No  Dither  (ME/I) 
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Figure  A-7.  7-Point  Parameter  Space:  Parameter  Variation,  No  Dither  (MIVA) 
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Figure  A-8.  7-Point  Parameter  Space:  Parameter  Variation,  No  Dither  (ME/I) 


Figure  A-9.  7-Point  Parameter  Space:  Swept  Sine  Wave  Dither,  29-33  Hz  (ME/A) 
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Figure  A- 10.  7-Point  Parameter  Space:  Sine  Wave  Dither  (ME/A) 
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Figure  A- 12.  7-Point  Parameter  Space:  Square  Wave  Dither  (MWA) 
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Figure  A-14.  7-Point  Parameter  Space:  Wide  Band  Noise  Dither,  7-33  Hz  (MBA) 


Figure  A- 15.  7-Point  Parameter  Space:  Wide  Band  Noise  Dither,  30-33  Hz  (MWA) 


Figure  A- 16.  7-Point  Parameter  Space:  Wide  Band  Noise  Dither,  30-33  Hz  (MB/I) 
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Figure  A- 18.  7-Point  Parameter  Space:  Rigid-body  Dither  (ME/I) 


Figure  A- 19.  29-Point  Parameter  Space:  No  Dither  (ME/A) 
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Figure  A- 20.  29-Point  Parameter  Space:  No  Dither  (ME/I) 


Figure  A-21.  29-Point  Parameter  Space:  Sine  Wave  Dither  (ME/A) 


Figure  A-22.  29-Point  Parameter  Space:  Sine  Wave  Dither  (ME/I) 
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Figure  A-23.  29-Point  Parameter  Space:  Square  Wave  Dither  (MHVA) 


Figure  A-24.  29-Point  Parameter  Space:  Square  Wave  Dither  (ME/I) 


Figure  A-25.  29-Point  Parameter  Space:  Wide  Band  Noise  Dither,  30-33  Hz  (MEJA) 
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Figure  A-26.  29-Point  Parameter  Space:  Wide  Band  Noise  Dither,  30-33  Hz  (MBA) 
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Figure  A-27.  29-Point  Parameter  Space:  Rigid-body  Dither  (ME/A) 


Figure  A-28.  29-Point  Parameter  Space:  Rigid-body  Dither  (ME/I) 
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